Glacier surface-area changes in Sagarmatha national park, Nepal, in the second half of the 20th century, by comparison of historical maps

Abstract

We investigate variations in the surface area of glaciers in Sagarmatha national park, Nepal, during the second half of the 20th century through comparison of a map applicable to the late 1950s with the official map of Nepal in the early 1990s. The comparison reveals a slight overall decrease in glacier area (by 4.9%, from 403.9 to 384.6 km²), a result which, though potentially subject to errors arising from cartographic interpretation, is in line with the area reductions found by other studies of Asian glaciers. We find that the areas of some individual glaciers, the largest situated at higher altitudes, increased during the study period. This was most apparent for the glaciers oriented to the south, with the increase occurring mainly in the glacier accumulation zones while the fronts tended to recede. Meanwhile, the smaller glaciers, situated lower and on steep basins, experienced a reduction. For the smaller glaciers, the sections most affected by change were the accumulation zones, and these glaciers showed a tendency for the front to advance. In this region there is a lack of climate data for high altitudes. Nevertheless, observations from stations situated around the park suggest that, alongside temperature variations which are often considered the primary factor eliciting glacier response, changes in precipitation play a significant role.

Introduction

Global climate change is manifested in mountain areas by a series of effects, magnified by the fragility which characterizes high-elevation ecosystems (Beniston, 2003). The great Asiatic mountain chains such as the Himalaya contain some of the largest ice masses outside the polar regions. These ice masses constitute a water resource which assures the survival of around 500 million people who inhabit the basins of the Indus, Ganges and Brahmaputra rivers (S. Sharma, http://www.water-2001.de). It is thus important to study remote mountain areas such as Sagarmatha national park (SNP), Nepal, which are subject to extreme climate conditions. Due to their biotic and abiotic characteristics, their ecosystems are highly sensitive and vulnerable and, for this reason, universally considered ideal sites for the study of long-term environmental change. Because variations in glaciers and lakes are strongly dependent on climate, they can be considered two very significant indicators of climate trends (Haeberli 1990; Wood 1990). Therefore, defining the current dimensions and recent variations of glaciers and lakes is important from the perspective of the water resources they represent and the climate changes they indicate. In the past few years, such measurements have been carried out in a number of studies using satellite images. Among these are the GLIMS project (Global Land Ice Measurements from Space) (Kargel and others, 2005) and some significant work concerning recent glacier variations in the high mountain ranges of Asia (Ye and others, 2006a, b; Berthier and others, 2007). However, these studies have only examined the evolution of glaciers since the 1970s, when satellite data became available.

The present work attempts to measure glacier area variations from an earlier date, the end of the 1950s, using a topographic map constructed at that time. The second map, utilized to achieve a comparison of glacier areas based on two data sources of the same type, is the latest available topographic map of SNP. These two maps allow us to observe the changes between the end of the 1950s and the beginning of the 1990s. The results of this analysis can be considered an initial step in tracking t`he evolution of glacier areas in SNP. In future, the comparison might be extended to the present-day configuration of the glaciers, subject to careful assessment of the comparability of the topographic maps used in this work with satellite images available for SNP (Bolch and Kamp, 2006; Buchroithner and others, 2006).

The present analysis of the evolution of water bodies in the Himalayan region is part of a project for the development of a Glaciological Information System, incorporated into a broader GEoDataBase of SNP, intended for collecting morphometric and climate data and tracking their evolution, spatially and temporally. The GeoDataBase for SNP is freely available for download from the website www.hkkhpartnership.org, thus making accessible a tool which comprehensively covers all the disciplinary fields relevant to climate change.

Study site

Sagarmatha national park is situated in the Solu–Khumbu district, in the northeastern region of Nepal (Fig. 1), and occupies the northernmost part of the Dudh Koshi river basin, which in turn is part of the Koshi river basin (or Sapta Koshi basin), one of the seven major hydrographic basins into which Nepal is subdivided. SNP, officially instituted on 19 July 1976, covers an area of 1141 km² and has unique geographical features, being surrounded on all sides by the highest mountain ranges on Earth. The terrain is extremely irregular, with altitudes ranging from 2845 m (Monjo) to 8848 m (Sagarmatha is the Nepalese name for Qomolangma (Mount Everest)).

Fig. 1. Location of the historical meteorological stations near Sagarmatha national park (SNP).

Byers (2005) describes the climatic features of SNP: geographically, it lies within the subtropical Asian monsoon zone characterized by pronounced summer rainfall maxima, with over 80% of the annual precipitation falling during an approximately 4 month period between June and September (Mani, 1981; Barry and Chorley, 1982). Winters are normally dry, although occasional mid-latitude cyclones, driven by the subtropical jet stream that takes its winter position just south of the Himalayan ridge, can cause heavy snowfall events (Klaus 1966; Zimmermann and others, 1986). A gradual transition from the dry winter season to the summer monsoon results from pre-monsoon convective precipitation events which are often accompanied by thunderstorms (Ueno and others, 1993). The prevailing axis of the monsoons is southern and southwestern (Müller, 1980). Ageta (1976) observes that while monsoon activity from the southern foot of the Himalaya decreases, approaching the interior of the main range of the Himalaya precipitation increases at some of the higher points. For instance, the average annual precipitation of Namche Bazar (27.83° N, 86.72° E; 3450 m) is ∼1000 mm a⁻¹ and decreases with elevation (to ∼500 mm a⁻¹ in Dingboche (27.89° N, 86.83° E; 4355 m), while the total amount of precipitation around peaks and ridges can be four or five times greater than that in valleys (Higuchi and others, 1982). Precipitation occurs mainly in the daytime on the ridges, but at night in the valleys.

The glaciers of SNP, like those of the Himalaya in general, exhibit certain morphological features of particular interest. Nearly all (28 out of a total of 29 in SNP) are ‘black glaciers’, known also as D-type or debris-covered glaciers; these are glaciers in which the ablation zone is almost entirely covered by surface debris that significantly alters the energy exchanges between the ice and the atmosphere. According to Moribayashi (1974), Inoue (1977) and Smiraglia (1998), these Himalayan glaciers are further characterized by a marked difference between their maximum elevation and the snowline elevation, a stable correlation between the accumulation area and the overall length of the glacier, abundant accumulation driven by wind and avalanches and multiple confluences that sometimes make it difficult to identify the principal lobe. Another characteristic feature is the debris cover which gradually increases toward the lower part of the glacier’s ablation zone. Here, the coverage becomes complete, and the glacier is regarded as inactive, presenting an apparently chaotic arrangement of detritus, depressions and cavities that often contain very small lakes of muddy water up to about 10 m in depth (Iwata and others, 1980). Fujii and Higuchi (1977) report that the C-type glaciers (i.e. glaciers without debris cover on their ablation zone) of SNP at the end of the 1950s occupied an average area of 0.43 km² (standard deviation 0.55 km²). Because the minimum threshold for considering glacier areas in this study was set at 1 km², as we shall discuss below, the only debris-free glacier we consider is Langmuche glacier.

Data and Methods

Cartographic supports and georeferencing

The topographic and thematic maps available for Nepal are well documented by T.B. Pradhananga (http://www.gisdevelopment.net/aars/acrs/2002/mgt), who provides a complete report of all the existing (topographical, land resource, thematic, derived, etc.) maps by Nepal’s Topographical Survey Branch. Of the existing cartography, only the Khumbu Himal map, published in 1978 and the current official map of Nepal published in 1997 are sufficiently representative of the area of SNP. Table 1 details the characteristics of these two maps, while Table 2 gives the nomenclature of the glaciers.

Table 1. Map features used for the analysis of glacier area variations

Table 2. Nomenclature of glaciers according to the different sources

Both maps have the same scale, so it is possible to compare the level of detail of the information they contain. Compared to the Khumbu Himal map, the official map of Nepal provides slightly less definition of qualitative aspects such as the morphological representation of rocks and, of more importance for this study, no information on debris cover. This latter shortcoming means we cannot directly evaluate the variations over time in the areas of the accumulation zones, generally without debris cover, and of the ablation zones, where debris is present.

The Khumbu Himal map covers 98.0% of the park’s territory, extending from longitude 86°33′ to 87°07′ E, while the boundaries of SNP extend in longitude from 86°31′ to 87°01′ E. However, this map does cover the entire latitude range of the park. To take this into account and permit comparisons with the current official map of Nepal, the data referred to the latter map have also been restricted, in the analysis and the tables that follow, to the longitudinal limit (86°33′ E) of the Khumbu Himal map.

In order to compare the glacier areas it is first necessary to define the time periods to which the maps refer. As Table 1 shows, the Khumbu Himal map takes in a wide span of time between the photographic survey (1921), the terrestrial photogrammetric survey (1935 and 1939) and the terrestrial survey, itself done over a lengthy period, 1955–63 (Schneider, 1967). To narrow down the dating of this map, we use the fact that Müller, as we describe below, attributes the data of his inventory conducted on the Khumbu Himal map (Schneider, 1967) and the subsequent analyses to ‘the 1950s and early 1960s’ (Müller,1980). On the strength of the period given by Müller, we consider the Khumbu Himal map to be representative of the late 1950s (and hereinafter refer to it as the 1950s map). The period of the current official map of Nepal is better defined, and can be considered representative of the beginning of the 1990s (and is hereinafter referred to as the 1990s map). A comparison of the two maps thus allows us to analyze the changes occurring over a period of nearly 35 years.

The two maps were georeferenced according to the projection system of the official map of Nepal. However this was not done by simply altering the projection system of the Khumbu Himal map. In fact, J. Gspurning and others (http://avalanchemapping.org/Repapers/GISinhighmountainKumbu.pdf) have pointed out the horizontal deviations which exist between identical map control points, due to the poor accuracy of the local projection systems. We therefore proceeded to resample the Khumbu Himal map using the ArcView^© tool, based on >100 map control points (the mountain peaks) evenly distributed on the maps. The a posteriori check showed that the maximum horizontal deviation (along the x and y axes) did not exceed 15 m (georeferencing error (GE)), compared with 336 m for the y axis and 19 m for the x axis observed by Gspurning and others.

Evaluation of the accuracy of the maps

Before analyzing the precision errors arising from the process of cartographic representation and its subsequent conversion from analogue to digital form, we must consider the problem of the interpretation of glacial features (accuracy). This factor must be taken into account to ensure that any differences observed between the two maps effectively reflect real changes, rather than artefacts of different degrees of accuracy in the construction of the two maps. For the specific purposes of this work, this implies evaluating the comparability of the glacier outlines on the two maps. This is unfortunately not a simple matter, as it would require retracing the process by which the cartographers discriminated between what is and what is not mapped as a glacial feature. Firstly, there are a multiplicity of data sources to consider: photographic material, photogrammetric images, direct observations, as well as the personal experience of the cartographer in interpreting the features on the terrain. Secondly, any assessment would have to include various forms of ice cover in order to be representative of the entire park area. Finally, we have to consider that the photogram-metric images used for constructing the Khumbu Himal map were taken on glass plates whose current state of preservation does not allow them to be easily re-examined.

In metrology, the interpretation of glacial features on maps is commonly evaluated in terms of a data-accuracy index, defined as the discrepancy between the measured value (mapped glacier outlines) and the real value of the studied quantity (actual glacier outlines). Such a discrepancy generally arises from biases introduced by the instrumental method or, as in our case, by human factors (interpretation of glacier features). The error is considered systematic because it remains constant on repeating the measurement, making it difficult to determine not only its magnitude but even its very existence, except through a repetition of the entire measurement process (which is not feasible in our case) or a consistency evaluation of the results (Glosup and Axelrod, 1998). Consequently, despite the difficulty of directly analyzing how the glacial features were interpreted on the two maps, it is important to determine whether the results of the comparison can be ascribed to real changes, so in the discussion of the results we attempt to evaluate the congruence between the results and the processes that may have brought them about.

Calculation of the cartographic error

A consideration of fundamental importance for Geographical Information System (GIS) processing is the precision of data. In the present work, the data processing involves conversion of topographic maps from analogue to digital form, entailing a series of approximations, with inevitable effects on data quality. A final cartographical representation is dependent on various factors, including the precision of the map, the georeferencing error and the precision of the operator, which gives rise to the vectorization error. Assuming the vectorization error to be minimal, and having dealt with the GE of the maps as discussed above (15 m along the x and y axes), we can now further account for the precision of the data.

Every map has a well-defined metric precision which derives from the process by which it was constructed. The precision of the location of points and the level of detail of a topographic map, whether in digital or analogue format, are determined by its scale factor. The precision of a derived cartographic representation is the degree of planimetric and altimetric tolerance between the coordinates of any given point on the source map and the coordinates of that same point on the derived map. Evaluating the precision of a derived map is therefore necessarily dependent on an evaluation of the precision of the source map. The graphical resolution limit of a map is assigned an arbitrary value of ∼0.2 mm, based on the threshold visible to the human eye. The level of approximation of a map is conventionally given by this limit of 0.2 mm multiplied by the scale factor. Therefore, for the source maps used in this work (scale 1 : 50 000), the approximation, i.e. the linear resolution error (LRE), corresponds to 0.2 mm × 50 000 = 10 m (Inghilleri, 1974).

The overall linear error (LE) associated with the distance between two points must take into account the GE discussed earlier, as well as the LRE, and is the root mean square (rms) of the two errors (Bevington and Robinson, 1992):

(1)

GIS-based processing is often used for analyzing variations in areas and volumes. It is necessary to determine the deviation between the actual values and those computed based on the digitized data. Because the perimeter of an area is delineated by manual digitizing, the value of the planimetric vectorization error is often used to quantify the areal resolution error (ARE), which indicates the range of tolerance of the area measurements in question. The ARE is usually computed as the product of the perimeter, l, of the area in question and the LRE (Inghilleri, 1974). Since the LE is greater than the LRE, the areal error (AE) associated with an area must therefore be computed as:

(2)

In Table 3, the AE is referred to the glacier area variations between the end of the 1950s and the early 1990s, computed as the rms of the AEs associated with the 1950s and 1990s areas. The overall error was also referred to the 1950s area of each glacier, thus directly compared with the glacier area difference.

Table 3. Glacier comparison table

The digital elevation model

The digital elevation model (DEM) was constructed by digitizing the contours of the 1990s map. The interpolation was done by the kriging method using the Surfer^© software, setting the north–south isotropy of the data identified from analysis of the semivariogram. A DEM pixel dimension of 20 m × 20 m was used. Considering that the mean elevation range of the glaciers is ∼1620 m for a mean length of ∼6.7 km, there were on average ∼8 cells interpolated between adjacent contour lines. To judge the validity of the above processing method and hence the precision of the morphometric data thus obtained, we computed a posteriori the average height deviation (along the z axis) between the interpolated values and ∼100 map control points, obtaining an average absolute difference of 6.3 m. The mean error associated with the average glacier elevation (Table 2) is therefore 0.2%. For the surface gradient, because the calculation is based on the relative values between consecutive DEM pixels, the accuracy of the DEM absolute values is less important (Quincey and others, 2007) and has here been disregarded.

Statistical analysis

The statistical method used in this work is a correlation analysis whose results are expressed as a coefficient of determination, R ². The frequency distributions of the dependent variables are also given to show that all correlations are based on Gaussian distributions. The correlations found are validated by a parametric inferential analysis of the observed correlation coefficient, using a Student’s t test (Venables and Ripley, 2002). This test evaluates the probability of the null hypothesis H₀: r = 0 (i.e. absence of linear correlation between the two examined variables, meaning that the correlations found are due to chance), and hence the probability that the correlations found are significant (H₁: r ≠ 0). A significance level (α) is generally chosen to define the rejection range of the null hypothesis as α = 0.05, α = 0. 01 or α = 0.001. The p value gives a measure of the plausibility of the null hypothesis, defined as the minimum significance level, α, of the test for which the null hypothesis would be rejected. If p < α the null hypothesis is rejected and the correlation is deemed significant. The significance of the correlation coefficient, r, is obtained using the table of critical values of r_{α ,n} , as a function of the number of observations, n, and the different significance levels, α. For the linear correlations discussed below, the observed correlation coefficient, r, was compared with its critical values, r_{α ,n.} If r> r_{α ,n} then p < α, and the null hypothesis was rejected. All the graphs included in this work were found to be statistically significant, and in each case we give the value of the coefficient of determination, R ², and the significance level of the associated correlation coefficient (p < α).

The meteorological data

There are not many weather stations above 2000 m in Nepal, as noted by Shrestha and others (1999). They provide a list of stations operated by the Department of Hydrology and Metereology (DHM), some of which have been active since the mid-1960s. Temperature and precipitation data are available as monthly means (Nepal: DHM, 1998). Figure 1 shows the locations of stations managed by the DHM situated near SNP. Also included in the figure is Shegar station (28.63° N, 87.08° E; 4300 m a.s.l.) in Tibet, for which annual temperature and precipitation graphs are described by Jin and others (2005). In the discussion we present annual mean temperature and precipitation graphs for Kathmandu station, as this recorder is the meteorological reference for Nepal. In this case, after a 9 year period of joint operation, a station located at the Indian Embassy (IE) (27.44° N, 85.20°E; 1324 m a.s.l.) was replaced by one located at Kathmandu International Airport (IA) (27.42° N, 85.22° E; 1336 m a.s.l.). We developed a linear regression model of the annual mean temperature and annual total precipitation data for the IE station, based on the IA annual data for the overlap period (1968–75). The extension of the IE record using the regression model (Kathmandu-simulated series) allows us to estimate the annual temperature and precipitation trends during the observation period. The error arising from the estimates is the rms of the differences between the reconstructed values and the observed values during the overlap period (Hays and Winkler, 1970).

Results

Glacier identification

Previous work has constructed inventories of the glaciers in SNP. Müller (1970) reports an inventory of glaciers in the Qomolangma region (corresponding in extent to the current SNP) which can be considered one of the pilot studies in this field. This work was undertaken using the map of the late 1950s on a 1 : 50 000 scale (Schneider, 1967). Fujii and Higuchi (1977), referring to this inventory, analyzed the relations between glacier formations and the topographic conditions of the glacier basins. Unfortunately for the purposes of the present study, it was not possible to use the area data from the above inventory directly, because it restricts the area calculation to those glaciers classified as ‘extensive’, giving only the geographical position for the smaller glaciers. Such a method results in a substantial underestimation of the total SNP glacier area (287 km²). In 2001 the International Centre for Integrated Mountain Development (ICIMOD; Mool and others, 2001) constructed a glacier inventory for all Nepal, and thus also for SNP. Working at the national level, a variety of data sources were used including satellite images, aerial photographs and maps. The total SNP glacier area was calculated to be 360 km². However, this inventory also proved unsuitable as a basis of area comparison for the present work, because the multiplicity of its sources cannot be attributed to a precise historical period. That said, since it is not the purpose of this work to construct a new inventory, but only to evaluate changes in glacier areas, it was decided to adopt the glacier codes proposed by ICIMOD and shown in Table 2. We see in Table 2 that glacier nomenclature differs considerably across the two maps. We choose to use the ICIMOD inventory as a reference, giving it preference in assigning unambiguous place names. For those glaciers that were not assigned a name in the ICIMOD inventory, we used, if reported, the name given in the maps (e.g. Phunki, Machhermo, Kyajo, Chhuitingpo) or the ICIMOD glacier code (e.g. Kdu_gr 181, Kdu_gr 125, Kdu_gr 38). The only remaining problem with the ICIMOD inventory concerned Ama Dablam glacier, which was swapped with Ombigaichain glacier, and has here been called Duwo glacier in accordance with both maps. Table 2 also shows the coordinates of the glacier centroids as identified on the 1990s map. Figure 2 shows the glacier areas at the end of the 1950s and in the early 1990s.

Fig. 2. Glacier comparison map of SNP.

For classification criteria used to define the extent of each glacier, we referred to the glacier classification guidance for the GLIMS glacier inventory (Rau and others, 2005). This gives different classification methods, depending on the purpose. For analyzing the extent of glaciers, the ‘Form parameter group’, which essentially describes the outline of the glacier, is suggested. Following the criteria given for this category, we associated glaciers with their catchment areas, thereby defining compound basins (fig. 20 of Rau and others, 2005) in which several individual accumulation basins feed into a single glacier system. In accordance with Mayer and others (2006), the tributary valley glaciers are included as part of the main glacier even if they are not directly connected to the accumulation and ablation zone of the main body in one of the two maps. Proceeding in this way, these areas were considered to be the glacier’s possible ‘feeding’ zone. Using the above classification criteria also ensures that the comparison is being made between the same glacier surfaces, without mistaking the joining or detaching of tributary glaciers as increases or decreases in glacier area. For instance, on the basis of these criteria, Gyubanare glacier, which is usually considered independent (Müller, 1970; Fujii and Higuchi, 1977), is not detached from Ngojumba glacier.

Working on the early-1990s map, we began by digitizing the surface hydrography of SNP and subsequently the corresponding hydrographic basins. We then proceeded to digitize the glacier areas on both the 1950s and 1990s maps. As discussed above, to permit a comparison between individual glaciers on the two maps, it was necessary to identify glacier basins, choosing as the basin outlet the glacier front which had advanced the most between the two periods under study. By so doing it was possible to ensure that each basin entirely contained the area of the corresponding glacier in both maps. To make comparisons between individual glaciers, we chose to consider only those glaciers whose area in the late-1950s map exceeded a threshold of 1 km². Proceeding in this way, a total of 29 glaciers, accounting for nearly the totality of glacier area in SNP (98%), were analyzed.

In Table 3 we see that the largest glacier was Ngojumba glacier, with a 1950s area of 98.7 km². In the 1950s the average elevation of the glacier areas was 5488 m, with Chhule glacier having the lowest (5099 m) and Khumbu glacier the highest average elevation (6100 m). The average aspect of the glaciers was 172°; in other words the glaciers were, on average, oriented southwards. Specifically, Tingbo and Cholotse glaciers were exposed to the west, Cholo glacier to the east and Kdu_gr 38 glacier to the northeast, while the remaining majority of glacier surfaces were instead exposed to the southwest. The average glacier slope during the 1950s was 28.0%. The two extremes were Kyajo glacier, with the lowest average slope (18.8%), and Phunki glacier, with the highest average slope (41.7%).

Table 3 also shows the morphometric characteristics of the glaciers in the 1990s. We observe that the total glacier area in the 1990s (without the 1 km² cut-off threshold) was 384.6 km² (AE = 20.7 km²), and the total glacier area beneath the threshold was 15.7 km²; the average glacier (above the threshold) area was 13.5 km², but with standard deviation (20.0 km²) nearly twice the mean. There is considerable variability in the size of individual glaciers, which can differ by nearly two orders of magnitude. In the 1990s, Ngojumba glacier was still the largest, with an area of 104.7 km². The average glacier surface elevation was 5491 m, with Cholo glacier having the lowest (5002 m) and Khumbu glacier the highest average elevation (6154 m). The average aspect computed in the 1990s was 176°. There were eight glaciers exposed to the south, six to the west, five to the southeast and five to the southwest. Langmuche, Langdak and Kdu_gr 38 glaciers were exposed to the northeast and Cholo and Thyangbo to the east. The average slope of the glaciers in the 1990s was 24.8%, with the two extremes being Kyajo glacier, with the lowest average slope (10.5%), and Phunki glacier, with the highest average slope (41.2%).

Glacier area variations

Figure 2 shows a visual comparison of glacier areas, while Table 3 gives the percentage change in the area of each glacier from the 1950s to the 1990s, and the morphometric parameters for these two periods. We observe that:

Overall, SNP experienced a net glacier cover reduction of 19.6 km² (4.9%) from 403.9 km² at the end of the 1950s to 384.6 km² at the start of the 1990s.

There were both positive and negative glacier area variations over the study period: in the 1990s, 22 out of 29 glaciers had decreased in area since the end of the 1950s (average loss of 26.4% per glacier, for a total of 32.0 km²), while 7 glaciers had increased in area (average gain of 10.2% per glacier, for a total of 17.0 km²). The greatest area increase was observed at Tingbo glacier, which expanded by 74.6%, while the greatest decrease was observed at Kdu_gr 38 glacier, which contracted by 53.5%. A visual analysis using Google Earth (http://earth.google.com) was conducted to identify any gross cartographic errors. This revealed that the enormous expansion of Tingbo glacier has to be considered an overestimate arising from an interpretation error in the 1990s map. Therefore the area changes of Tingbo glacier were excluded from further analyses.

The average slope of the glaciers decreased from 28.0% to 24.8% over the study period. The average aspect of the glaciers changed from 172° to 176°.

In summary, our comparisons revealed that during the study period SNP experienced a small net reduction in glacier cover. The changes to individual glacier areas were variable, with a reduction in size of a large number of glaciers being offset by an increase in size of a smaller number of glaciers, thus mitigating the overall loss of glacier area. In the following subsection, we examine these opposing variations in more depth, attempting to correlate them with the morphological characteristics obtained from the DEM.

For each glacier, the change in surface area from the 1990s to the 1950s is compared with its size (area) in the 1950s. Figure 3a shows the percentage change in glacier area as a function of glacier area. For the smaller glaciers (<20 km²) there is a clear tendency toward loss of glacier area over the study period, but no correlation between glacier size and area change is evident. The larger glaciers (>20 km²), however, are less likely to have lost area over the study period, and overall there is no clear trend toward loss or gain. In Figure 3b we see that the average area loss of the glaciers was 18.2%, although the effect of the larger-sized glaciers, which on average increased in area, meant that overall the park lost only 4.9% of its glacier area.

Fig. 3. (a) Glacier area variation from the 1950s to the 1990s, as a function of glacier size. (b) Frequency distribution of glacier area variation.

Figure 4a shows the percentage change in glacier area over the study period, as a function of glacier average elevation. Lower-altitude glaciers lost the most area, while higher-altitude glaciers increased their area (R ² = 0.34; p < 0.001). The glaciers of SNP had an average elevation of 5488 m in the 1950s, ranging from 5099 m (Chhule glacier) to 6100 m (Khumbu glacier) (Fig. 4b).

Fig. 4. (a) Glacier area variation from the 1950s to the 1990s, as a function of glacier elevation. (b) Frequency distribution of glacier mean elevations.

From Figure 5a we see an inverse relation (R ² = 0.21; p < 0.01) between percentage change in glacier area and average glacier slope. The glaciers which lost most area were those of highest slope, while the glaciers whose area increased were those of lower slope. The average slope of the glaciers was 28.6% in the 1950s, and the distribution is close to normal, but with three glaciers with slope >40% (Kdu_gr 38, Phunki and Tingbo glaciers) (Fig. 5b).

Fig. 5. (a) Glacier area variation from the 1950s to the 1990s, as a function of glacier slope. (b) Frequency distribution of glacier average slopes.

We also analyzed the relationship between glacier area and slope. From this comparison, it emerged that although there is, as expected, an inverse relation, its magnitude is not significant due to the presence in the park of several small glaciers with low slope gradients.

Figure 6a shows the glacier area in the 1950s relative to glacier aspect. The plot shows the larger glaciers were prevalently oriented to the south. The orientation of glaciers whose area decreased (empty circles) and increased (black squares) between the 1950s and the 1990s is indicated. Glaciers that increased are prevalently oriented toward the south. Figure 6b and Table 4 highlight this distribution, showing the area variation of each glacier (absolute value in km²). The glaciers which expanded were oriented to the south, while those which contracted were oriented toward other sectors. The greatest reductions were observed in the southwest sector, and amounted to 10.4 km², representing 39.0% of the total glacier area losses.

Fig. 6. (a) Glacier size as a function of glacier aspect. The empty circles represent glaciers that experienced area loss between the 1950s and the 1990s, while the black squares represent the glaciers whose surface area increased. (b) Grouping of total glacier area variations (km², radial axis) according to their orientation sector (°, circular axis).

Table 4. Grouping of total glacier area variations according to glacier aspects. The total differences are plotted in Figure 6b on the radial axis. They represent the area changes experienced between the 1950s and the 1990s by all the glaciers grouped within the same orientation sector, while the negative and positive variations are obtained by considering exclusively the glaciers which have decreased or increased their area

Glacier area variations and changes in morphological features

Figure 7a shows the relationship between change in glacier area and change in glacier altitude during the study period. The glaciers which lost most area were, in the 1990s, situated at a lower average altitude than they were in the 1950s, while the glaciers which increased in area also exhibited an increase in average altitude (R ² = 0.29; p < 0.01). The area variations chiefly affected the accumulation zones of the glaciers rather than the ablation zones; in fact if the expanding glaciers increased their average altitude, they must necessarily have expanded more in their accumulation zones than in their ablation zones. Conversely, the contracting glaciers, whose average altitude decreased, must necessarily have suffered greater losses in their accumulation zones than in their ablation zones. We note, in any case, that the average elevation change (−1 m) is very small (Fig. 7b).

Fig. 7. (a) Glacier area variation from the 1950s to the 1990s, as a function of change in average elevation. (b) Frequency distribution of glacier elevation differences between the 1990s and 1950s.

Very similar behaviour is shown in Figure 8a, a plot of glacier area change over the study period against glacier slope change. Glaciers which shrank the most also showed a reduction in slope gradient, while the glaciers that expanded exhibited increased slope (R ² = 0.46; p < 0.001). The mean change was a reduction in slope of 0.034 (Fig. 8b).

Fig. 8. (a) Glacier area variation from the 1950s to the 1990s, as a function of change in slope. (b) Frequency distribution of glacier slope differences between the 1990s and 1950s.

These last two findings support the hypothesis that glacier area change, whether positive or negative, is primarily the result of changes occurring in the accumulation zones rather than in the ablation zones. This allows us to overcome the limitations of the cartographic representations. (As noted previously, the official map of Nepal, unlike the Khumbu Himal map, does not record any distinction between ablation zones (generally debris-covered) and accumulation zones (generally free of debris).)

The hypothesis that the greatest glacier area changes occur in the accumulation zones does not rule out the possibility that notable variations may also occur at the glacier front. Figure 9b is a map of Ngojumba glacier, the largest glacier in SNP, which expanded from 98.7 km² in the 1950s to 104.7 km² in the 1990s (+6.1 %, AE = 3.4%). Figure 2 shows its location in the north of SNP and its south-facing aspect (182°). We see that the area of the accumulation zone has increased. Figure 9a is a map of Cholo glacier, one of the smallest in SNP, which shrank from 2.1 km² in the 1950s to 1.5 km² in the 1990s (−31.2%, AE = 10.1 %). Figure 2 shows its location in the centre of the park and its east-facing aspect (83°). More specifically, Cholo glacier is situated in a western valley of the Khumbu valley and drains into the Imja Khola river. We observe here that the area of the accumulation zone has substantially decreased.

Fig. 9. Two examples of glaciers in which the greatest positive and negative area variations occurred primarily in the accumulation zones: (a) Cholo glacier, (b) Ngojumba glacier.

Discussion

The above analysis allows us to compare the glacier area changes in SNP on the basis of two ‘snapshots’, at the end of the 1950s and in the early 1990s. The overall reduction in glacier area within SNP was found to be 4.9%, but with a margin of error associated with the cartographic interpretation (AE) of 4.9%. This uncertainty, according to Mi and others (2002), is typical of glacier area determinations based on topographic maps (∼5%). Although with a wide margin of error, other studies, integrating topographic maps and remote-sensing data, have also found reductions of the same magnitude in some Asian regions. Table 5 gives details of studies conducted on glacier area changes over a comparable timespan and covering an extensive area. Although these regions are very likely to be subject to different climate conditions and include glaciers of different types, the overall reduction in glacier area reported by the various authors is in line with that found in the present work. The consistency with these other results indicates that, even though it was not possible to re-analyze the processes leading to the interpretation and cartographic representation of the glacier outlines, the observed variations in glacier area are plausible.

Table 5. Variations in glacier areas in other Asian regions during a comparable historical period to that analyzed in the present work. TM: Thematic Mapper; MSS: multispectral scanner

Liu and others (2002) performed a comparative analysis of glacier area change in the northern part of the West Qilian Mountains (northwest China) between 1956 and 1990. Their results indicate a 4.8% reduction in glacier area. Ye and others (2006a, b), studying the western Himalaya and central Tibet, found reductions in glacier area of 2.8% (between 1976 and 1990) and 2.3% (between 1969 and 1992), respectively. Karma and others (2003) studied 66 C-type glaciers in the Bhutan Himalaya by comparing topographic maps from 1963 with satellite images from 1993 and found that the glaciers had retreated by 8.1%. These studies additionally found an acceleration in glacier retreat, starting in the 1990s.

Glaciers change as a consequence of climate fluctuations. Among many factors, summer temperature and snowfall are the two most important variables controlling glacier development on the Himalaya. Summer temperature dictates ablation, while snowfall influences accumulation. Investigations have shown that the increase in air temperature during the last century has started to affect glaciers in the Himalaya. Shrestha and others (1999) studied relationships between temperature trends in Nepal and global-scale trends by comparing the Kathmandu annual mean temperature with the 24–40° N annual temperature anomaly. They found some striking similarities between these records. The 24–40° N annual temperature shows a decreasing trend from 1938 to 1972, and an increasing trend from 1972 to 2000. In the Kathmandu record there is a general cooling trend from 1935 to 1974, followed by subsequent warming, but with different magnitudes than the 24–40° N record. The cooling from 1938 to 1970 lowered the 24–40° N annual mean temperature by ∼0.23°C, while in the Kathmandu record the drop in annual mean temperature was ∼0.6°C. Similarly, the warming after 1970 raised the 24–40° N annual mean temperature by ∼.5°C, while for Kathmandu the annual mean temperature rose by ∼1°C over the ∼30 year period. Figure 10 shows the annual mean temperature trends for the stations situated near SNP, and the Kathmandu-simulated series from the end of the 1950s to the early 1990s. All the stations show an increasing temperature trend over the full study period, even though there is no complete dataset available for the series closest to SNP. Table 6 shows, for each station, the mean annual temperature and the mean increases per year over the periods for which data are available. The same table also illustrates the annual mean increase taking into account only the temperatures recorded during the summer season (June–September), which, as previously noted, are one of the most important variables controlling glacier development in the Himalaya. We see that, in some cases, greater temperature increases were recorded during this period.

Fig. 10. Trend of annual mean temperature for the stations situated near SNP and for the Kathmandu series obtained as described in the text. The Shegar station graph was produced by Jin and others (2005). Data for the other stations are from DHM (1998) and were processed for this work.

Table 6. Annual mean temperatures, mean increase per year and mean increase per year considering exclusively the monsoon season (June–September) observed during the period for which data are available for each station

It is possible that the period of global cooling between 1940 and 1970 may have resulted in a glacier readvance. Although the effects of gradual global warming during the 20th century have been clearly manifested in the cryosphere, it is also true that this has not been a linear process, with numerous slowdowns and even reversals (Wood, 1988). This has also been recorded between the 1960s and 1980 in some glaciers of the Karakoram (Smiraglia, 1998). Goudie and others (1984) and Gardner (1986) observed a reduced rate of recession in this region and, in some cases, advances, noted for instance on Biafo and Rakhiot glaciers.

Taking into account the variability of the temperature trends during the study period, the delayed response of debris-covered glaciers to temperature variations and the tendency in debris-covered glaciers for thickness to decrease, rather than surface area (down-wasting) (Kadota and others, 2000), it is not surprising that overall glacier area is only slightly reduced. The D-type glaciers of the Himalaya are not capable of dynamically adjusting to an accelerated warming by retreat. Rather, they react by down-wasting and decoupling of glacier parts. Nakawo and others (1999) describe recent changes in the debris-covered area and the relevant thickness of Khumbu glacier. In many cases, down-wasting is associated with rapid growth in the number and dimensions of supraglacial lakes. These can generate catastrophic glacier lake outburst floods (GLOFs), which constitute the strongest manifestation of the deglaciation process and its attendant risks in the region.

Divergence between individual glacier area variations

The glaciers of larger size, which were situated at higher altitudes, increased in size during the study period. Conversely, the smaller glaciers, situated at lower altitudes, showed a reduction in size (Figs 3 and 4). Similar findings by Paul and others (2004) on a sample of 938 alpine glaciers from 1973 to 1998/99 and Jin and others (2005) for the Tibetan Plateau (∼900 glaciers) show a larger contribution by small glaciers to the total loss of area during the study period. The higher sensitivity of small glaciers to climate effects can be explained by considering that when the equilibrium-line altitude increases, it can extend above the altitude of many of the small glaciers, whereas the debris on the tongues of large glaciers often prevents them melting (Jin and others, 2005). This result indicates that the interpretation of the glacier outlines on the two maps can be considered realistic, as it is consistent with these other studies.

We observed that area increase was most apparent for the glaciers oriented toward the south (Fig. 6b) and occurred principally in the glacier accumulation zones (Figs 7 and 8). From Figures 2 and 5a we noted also that the larger glaciers (Ngojumba, +6.1%, AE 3.4%; Bothe Kosi, +12.6%, AE 3.4%; Lumsamba, +7.6%, AE 4.4%; and Nuptse, +25.4%, AE 7.3%), i.e. those exhibiting increased area, are predominantly oriented to the south. Conversely, area reductions during our study period were most evident in the glaciers oriented in other directions. Also in these cases, the area of the glaciers most affected was the accumulation zone. We need to determine whether the glacial features have been interpreted correctly on the two maps. Although there are a few studies (e.g. Hewitt, 2005) which report variations in accumulation zones that exceed those in ablation zones, these can plausibly be accounted for in terms of climate change, as we show below. The fact that the glaciers whose area decreased showed the largest changes in their accumulation zones is consistent with their steeper slopes and lack of debris cover. Melting induced by temperature rises may have caused more extensive decoupling, fragmentation and avalanches here than in the flatter, debris-covered ablation zones, which are more susceptible to down-wasting than to area loss.

The fact that the larger glaciers are primarily oriented to the south (Fig. 6a) is accounted for by the southern orientation of the main valleys in the area, and the predominantly south to southwesterly origin of precipitation patterns (Müller, 1980). The monsoon travels south–north in this area, before veering to the west (Barry, 1981). In addition, these glaciers are not occluded to the south by mountain ranges perpendicular to the monsoon that might act as a barrier against the mass of moist air which produces precipitation. If these glaciers have attained their present large size thanks to a favourable orientation with respect to monsoon precipitation, we conjecture that their increase in area may likewise have been favoured by an increase in precipitation. However, there was no increase in the area of glaciers situated at lower latitudes which, in addition to being smaller, have an orientation that diverges from the south–north axis. As shown in Figure 2, there are many glaciers which exhibit these characteristics (Chhule glacier, Nare glacier, etc.).

The above observations suggest that alongside temperature variations which are often considered to be the main cause of glacier change, variations in precipitation may also play a determining role. Temperature variations alone cannot account for why some glaciers increased in size while others (all of the same type) decreased in size. This latter observation is better explained by considering that there was also a precipitation increase over the study period. Figure 11 shows the annual precipitation trends for the stations situated near SNP and the Kathmandu-simulated series from the end of the 1950s to the early 1990s.

Fig. 11. Trend of total annual precipitation for the stations situated near SNP and for the Kathmandu series obtained as described in the text. For Shegar station, the graph was produced by Jin and others (2005). Data for the other stations are from Nepal: DHM (1998) and were processed for this work.

All the stations show an upward trend in precipitation, though for the series closest to SNP there is no complete dataset available. Table 7 shows, for each station, the yearly total precipitation for the periods for which data are available, and the corresponding observed increase per year over the period. It is not possible to make a direct comparison between these data, which come from different observation periods. We observe from Figure 11, however, that all the stations recorded a more pronounced increase until the mid-1970s, and a less steep rise in the following period. Here again, however, there remains uncertainty arising from the absence of data for the higher altitudes.

Table 7. Annual total precipitation and mean increase per year observed during the period for which data are available for each station

Glacier orientation relative to the prevailing direction of the monsoons may cause some glaciers to receive more precipitation than others, favouring those exposed to the south. What is more, increased precipitation may explain the increase in the area of the glacier accumulation zones. For these glaciers, the effects of the precipitation may predominate over effects of temperature, and the high elevations of these accumulation zones protect them from melting.

For the glaciers oriented in other directions and at lower elevations, which experienced a reduction in accumulation-zone area over our study period, the effects of temperature will have prevailed. Fukui and others (2007) show, testifying to the effect of temperature at these altitudes, that in the Khumbu glacier valley the permafrost lower limit has risen 100–300 m between 1973 and 1991, stabilizing at 5400–5500 m. These glaciers, in addition, are often very steep, with notable discontinuities. As Müller (1980) reported, this is a common feature of many glaciers in SNP: the extremely steep and large headwalls and side-walls of many valley glaciers contribute large amounts of avalanche snow to the accumulation area, frequently causing very irregular surface contours and, thus, irregular firn lines. An examination of Cholo glacier in Figure 9a shows a very steep glacier headwall in the 1990s, with remnant ice around the ridge-line at the highest elevation, and a disconnected lower-elevation glacier tongue, presumably fed by avalanches. For these glaciers, as explained above, it is therefore plausible that the accumulation zones, with steep slopes, may have lost more snow and glacial cover than the flatter, debris-covered glacier tongues, which instead have experienced greater loss in thickness (down-wasting) than in surface area (Nakawo and others, 1999).

A similar interpretation to that proposed above has been suggested in work relating to another period, which investigated the expansion of the glaciers of the Karakoram (Diolaiuti and others, 2003; Hewitt, 2005). For these glaciers, situated at high altitude, an increase in winter precipitation was observed. Effects of the rise in winter temperatures resulting from global warming were found to be less critical than winter precipitation increase for these high-altitude glaciers (Archer and Fowler, 2004). The glaciers at intermediate altitude, according to Hewitt (2005), did not benefit from this expansion and instead continued their process of retreat.

Conclusions

We have provided a comprehensive picture of the changes in glacier surface area occurring between the end of the 1950s and the early 1990s on the glaciers of SNP. An overall reduction in glacier area of 4.9% was observed, consistent with other studies that have found reductions of about the same magnitude in other Asian regions. However the area changes of individual glaciers were found to be variable, and temperature variations, despite being the primary cause eliciting glacier response, cannot alone account for why some glaciers increased while others decreased.

We have seen that between the end of the 1950s and the early 1990s the stations situated near SNP, in Tibet and Kathmandu, recorded increases in both temperature and precipitation. We propose that, in the case of the glaciers which expanded, the effects of precipitation predominated over those of temperature, owing to their higher elevations which protected them from melting, and favourable orientation for capturing monsoon precipitation. The area losses were instead experienced by the smaller glaciers, situated at lower altitudes, and with orientation less aligned with the prevailing precipitation where the effects of temperature dominated.

We conclude this analysis by noting that the main objective of the study was to examine the glacier area variations in SNP starting from the end of the 1950s, to which end a comparison of historical maps was carried out. It is important to point out the elements of uncertainty inherent in this type of analysis which cannot be overcome with the current state of knowledge. However, these analyses can serve as a prompt for further investigations. Furthermore, this study highlights the need for additional research into the effects which the observed changes may have caused, and for extending the analysis to reach the present-day configuration of the glaciers.

It should be remarked, in this connection, that the majority of studies on glacier area changes, as noted by Hewitt (2007), suffer from the absence of high-altitude data which makes it problematic to extrapolate the changes observed at weather stations to the glaciers. That said, starting from the early 1990s, there are data available from the weather station operated by Ev–K²–CNR at 5100 m which can be used to determine the cause of any changes between then and the present-day configuration of the glaciers.

Another important factor is the uncertainty arising from the accuracy of the maps used. Despite the difficulty of directly analyzing how the glacial features were interpreted on the two maps, the changes observed in this study do plausibly reflect real changes, as they are consistent with the processes that may have brought them about. However, the above considerations apply to the SNP glaciers taken together as a group, and detailed evaluation on the individual-glacier scale is likely to reveal further problems connected with the interpretation process. Many of the glaciers in SNP are characterized by steep slopes. This makes it very difficult, when creating a map, to determine exactly which parts of the basin are glacier, which are left-over avalanche debris and which are disconnected snowpatches.

In the near future it is planned to extend the analysis to the present day, to include the historical period when the greatest rises in temperature have occurred. Due to the absence of a more up-to-date cartography, we will use satellite images to complete the picture of the evolution of the glaciers.

Also planned is a study, complementary to the present work, which will include all the lakes in the region. Our aim is to more clearly determine the magnitude of the changes affecting water resources in the SNP region during the second half of the 20th century.