We evaluated elevation changes at four sites on debris-covered ablation area of Khumbu Glacier, Nepal Himalaya, since 1978. In 2004, we carried out a ground survey by differential GPS in the upper- and lowermost areas of the ablation area. The amount of surface lowering was calculated by comparing digital elevation models (DEMs) with 30-m grid size, as generated from survey data corrected in 1978, 1995, and in the present study. Because we could not access the middle parts of the debris-covered area due to surface roughness, for this area we used an ASTER-DEM calibrated by the ground survey data. The amount of surface lowering during the period 1978–2004 was insignificant near the terminus. A remarkable acceleration of surface lowering was found in the middle part of the debris-covered ablation area, where the glacier surface is highly undulating. In the uppermost area, surface lowering has continued at a steady rate. Surface flow speeds have decreased since 1956, revealing that the recent decrease in ice flux from the upper accumulation area would have accelerated the rate of surface lowering of the debris-covered area of Khumbu Glacier during the period 1995–2004.
Introduction
Changes in mountain glaciers are considered a reliable indicator of climate change (IPCC, 2007). In the Himalayas, few studies have reported changes in glacier mass based on field surveys (e.g. Fujita et al., 1997, 2001; Wagnon et al., 2007), whereas several studies have examined temporal changes in glacier length and area (e.g. Mayewski and Jeschke, 1979; Yamada et al., 1992; Kargel et al., 2005; Salerno et al., 2008). The area of glaciers covered with supraglacial debris accounts for about 36% of the total glacierized area in the Khumbu region of east Nepal (Fujii and Higuchi, 1977). The results of field experiments and observations reveal that ice under a thick mantle of debris melts more slowly than does bare ice, due to the insulation effect of the debris (e.g. Mattson et al., 1993). However, Sakai et al. (2000, 2002) showed that debris-covered glaciers experience considerable melting at exposed ice cliffs and where covered by supraglacial pond water, consistent with the ablation rate of debris-covered ice calculated from Landsat TM data by Nakawo and Rana (1999).
Because the existence of a debris mantle makes it difficult to detect temporal changes in the glacier area, it is useful to calculate changes in elevation when studying fluctuations in debris-covered glaciers. Digital elevation models (DEMs) obtained by photogrammetry or altimetry have made it possible to determine changes in the elevation of glacier surfaces worldwide (e.g. Arendt et al., 2002; Muskett et al., 2003; Rignot et al., 2003; Surazakov and Aizen, 2006). Although several recent studies have analyzed Himalayan glaciers based on remote sensing DEMs (Berthier et al., 2007; Bolch et al., 2008a), they did not perform on-site verification of the results. In contrast, changes in the elevation of parts of Khumbu Glacier have been surveyed by Watanabe et al. (1980), Kadota et al. (2000), and Kadota et al. (2002). These studies surveyed the ablation zone of Khumbu Glacier during the years 1978, 1995, and 1999, respectively, with Kadota et al. (2002) reporting accelerated surface lowering in the uppermost ablation zone during the period 1995–1999.
Changes in the elevation of a glacier surface reflect trends in ablation and emergence velocity (Paterson, 1994). In the ablation area, emergence velocity is defined as uplift of ice by convergence of ice flux which compensates lowering by ablation. The spatial distribution of surface flow velocities for glaciers in the Khumbu region has been estimated based on remote sensing data, including satellite radar interferometry, satellite radar feature-tracking, or image-matching between optical satellite images (Seko et al., 1998; Bolch et al., 2008b; Luckman et al., 2007; Quincey et al., 2009; Scherler et al., 2008). These studies have reported a reduction in ice flux in the middle part of Khumbu Glacier since the 1950s. However, little attention has been paid to error evaluation or bias calibration based on field survey data.
In this paper, we evaluate changes in the elevation of Khumbu Glacier since 1978, based on ground surveys during the period 1978–2004 and a DEM generated based on remote sensing data collected in 2004 by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). We also consider the mechanism that underlies the observed changes in elevation, based on a comparison between the contributions to glacier elevation of ablation and emergence velocity derived from flow velocity data.
Study Area
Khumbu Glacier flows down from the southwest face of Mt. Everest (8848 m a.s.l.) to its terminus at 4900 m a.s.l. along a length of about 16 km. The glacier consists of a debris-free accumulation area and a debris-covered ablation area, separated by an icefall located at an elevation of about 5600–5700 m a.s.l. The mean annual temperature of −3.0 °C and mean monsoon seasonal cumulative precipitation of 237 mm during the period 2006–2008 at Nepal Climate Observatory–Pyramid (5079 m a.s.l.), located on the right bank of Khumbu Glacier, were reported by Bonasoni et al. (2010). The majority of precipitation in the studied area occurs during the summer season (Ueno et al., 2001) under the influence of a monsoon. Hence, glacier accumulation and ablation in the Nepal Himalaya occur simultaneously during the summer monsoon season (Ageta and Higuchi, 1984). The contribution by avalanches to accumulation has been estimated to be double that by snowfall (Inoue, 1977). Watanabe et al. (1980) established four survey sites within the ablation area of the glacier (Areas 1 to 4; Fig. 1). Intermittent surveys at these sites were performed in 1995, 1999 (Kadota et al., 2000, 2002), and 2004 (this study). Area 1 is located near the terminus and is covered with debris more than 2 m thick, beneath which the glacier ice is considered to be stagnant (Nakawo et al., 1986). Areas 2 and 3 are characterized by rugged topography in the middle and lower parts of the ablation zone. Area 4 is located in the uppermost part of the debris-covered area, characterized by large ice pinnacles (Fig. 2).
Data
DIGITAL ELEVATION MODELS FROM GROUND SURVEY
Watanabe et al. (1980) generated contour maps of the ablation area in 1978 based on theodolite surveys at all four of the areas (Table 1). They included coordinate information for bench marks in Areas 1, 2, and 3. However, because Area 4 does not include such coordinate information, no contour map was compiled for Area 4 in 1978. Kadota et al. (2000) compiled contour maps of Areas 2, 3, and 4 for 1995 based on theodolite surveys, and constructed cross sections through Area 1 (Table 1). To generate DEMs from these earlier surveys, we digitized the 1978 maps of Areas 1–3, and the 1995 map of Areas 2 and 3. The contour interval for these maps is 5 m. In a similar way, we processed the 1995 contour map of Area 4 (10 m interval) and related point data (n = 35) measured for Area 1, as obtained for several survey transects (Watanabe et al., 1980; Kadota et al., 2000). Some of the earlier contour maps include spot elevations in addition to contour line data. DEMs were generated by regularized spline with tension (Mitášová and Mitáš, 1993) using open-source GRASS GIS software. The employed algorithm is able to process both line and point data in generating a DEM. The accuracies of the contour maps compiled using the method proposed by D'Agata and Zanutta (2007) varied from 1.66 to 3.34 m, as shown in Table 1 and as calculated as follows:
where σh is the root mean square error of DEM elevation, σγ is the map reading error (0.2 mm divided by the scale of the map), e is an empirical number commonly within the range 0.16–0.33, CI is the contour interval, and α is the local slope of the DEM. We used a conservative value of 0.33 for e to obtain the maximum error.
TABLE 1
Summary of survey data analyzed in the present study. The elevation accuracies of contour maps calculated based on the method proposed by D'Agata and Zanutta (2007).
In October 2004, we surveyed Areas 1 and 4 with a differential code-phase GPS (Nikon GPS finder pro-XR; accuracy 0.5 m) and theodolite with a laser distance meter (Nikon GF-205C). Individual surveys were co-registered by referring to benchmarks installed on bedrock near the glacier (Fig. 1). The network of benchmarks established in 1978 (Watanabe et al., 1980), has been maintained and extended to the present (Kadota et al., 2000). The relative positions of these benchmarks were surveyed by theodolite. We then co-registered the positions of the benchmarks (Universal Transverse Mercator zone 45N, WGS-84 datum).
The generated DEMs were set to a grid resolution of 30 m based on the results of test calculations using various resolutions (see below).
ASTER DEM
The ASTER sensor, which has stereo-pair capability in the Visible and Near-Infrared Radiometer (VNIR) band, is able to generate DEMs. Such ASTER-derived DEMs (ASTER DEMs) have been used previously to evaluate temporal variations in the dimensions of glaciers in the Everest region, Central Asia, and southern Patagonia (Bolch et al., 2008a; Khalsa et al., 2004; Rivera et al., 2005). We used a Level 4A1Z product distributed by the ASTER GDS at the Earth Remote Sensing Data Analysis Center (ERSDAC) in Japan. This semi-standard orthorectified DEM from the Level-1A data was produced using data from two telescopes: nadir-looking VNIR (band 3N) and backward-looking VNIR (band 3B). Two images of the same region, taken approximately 55 seconds apart, constitute a stereo pair. The ASTER DEM is automatically generated without ground control points by stereo photogrammetry. Details of the algorithm used for DEM generation can be found in Fujisada et al. (2005), Toutin (2008), and on the ERSDAC Web site (ERSDAC, 2002).
FLOW VELOCITY
Seko et al. (1998) evaluated surface flow velocities of Khumbu Glacier by tracking temporal changes in the locations of two gaps in the pinnacle zone in Area 4 (G1: upper gap, G2: lower gap; Fig. 3) using a 1956 map (Müller, 1959), a 1/5000 map compiled in 1978 by Iwata et al. (1980), a 1/50,000 map published in 1988 by National Geographic magazine (1988), and System Pour l'Observation de la Terre (SPOT) high-resolution visible (HRV) images collected in 1987 and 1995. We also calculated flow velocities between 1995 and 2004, as derived from survey data related to obvious boundaries between two ice pinnacles in Area 4. By averaging two points (the upper and lower sides of each pinnacle), we negated the influence of shrinkage of the pinnacles.
FIGURE 3
ASTER image (2004) of Area 4, showing the outline of ice pinnacles in 1995 (black line). Measurement points are shown at the upper (1) and lower (2) edges of pinnacles A and C (P-A1 and P-A2 for Pinnacle A; P-C1 and P-C2 for Pinnacle C). Two dashed orange lines (UECE and LECE) present the boundaries used for calculating a continuity equation. See the Method section for details.
Methods
DEM PROCESSING
A Level 4A1Z product of ASTER acquired on 10 November 2004 was used for Areas 2 and 3, because we could not conduct ground surveys in these areas due to problems with access. The Level 4A1Z product, which is automatically generated without ground control points, generally contains horizontal and vertical biases. We calibrated the horizontal biases using a gap-filled Shuttle Radar Topography Mission (SRTM) DEM (Jarvis et al., 2008) by minimizing the standard deviation of the elevation differences. The absolute horizontal accuracy for SRTM has been reported to be ±20 m, whereas the horizontal accuracy of DGPS measurements is 0.5 m. Hence, our coregistration of the ASTER DEM (adjusted to SRTM) and DGPS DEM is likely to have a bias of 1 grid as the worst case. This bias may cause a large erroneous uplift/lowering of the boundary between the side moraine and glacier surface, due to the steep slope at the boundary. However, The DEM differentiation of our analysis (Fig. 7) does not show such characteristics. Consequently, we consider that our results are not affected by the difference in horizontal accuracy between the ASTER DEM and DGPS DEM.
The elevations in the 2004 ASTER DEM were compared with those of the 2004 differential GPS-derived DEM (DGPS DEM) in Areas 1 and 4. Because the ASTER DEM and DGPS DEM data were acquired within 1 month of each other, changes in elevation between the DEMs related to ablation and emergence velocity are largely negligible. ASTER DEMs are known to contain errors regarding the elevation of narrow ridge-like features such as lateral moraines (Fujita et al., 2008). We performed DGPS measurements mainly on the glacier surface during the survey campaign in 2004. Hence, measurement data collected on valley floors, which represent a stable and planar topography, are insufficient for validation. However, given the relatively smooth nature of the glacier surface in Areas 1 and 4, we considered it appropriate to perform the validation using data corrected from the glacier surface. Therefore, we calculated the bias between ASTER DEM and DGPS DEM only for the glacier itself.
Figure 4 shows six elevation biases of the ASTER DEM compared with the DGPS DEM at three different resolutions (15, 30, and 90 m), originally provided by ASTER GDS. The result shows that the coarsest resolution (90 m) has the largest standard deviation. Changes in elevation of the glacier surface show the same trend at different resolutions. For subsequent analysis, we selected a grid resolution of 30 m, provided a balance between accuracy and calculation cost.
As shown in Figure 4, the elevation in the ASTER DEM is lower than that in the DGPS DEM by an average value of 14.5 m, with the difference being −14.0 m in Area 4 and −15.0 m in Area 1. A similar negative bias has also been reported for the Bhutan Himalaya (Fujita et al., 2008). We found that the standard deviation (SD) of vertical differences (7.2 m for Area 1, 7.1 m for Area 4; average value, 7.2 m) is smaller than the estimated accuracy of the ASTER DEM (estimated to be 10 m by Fujisada et al. [2005] and 11.0 m by Fujita et al. [2008]), even in areas of rugged topography. Hereafter, we set the bias and SD of vertical differences in the study area to values of −14.5 m and 7.2 m, respectively.
COMPONENTS OF SURFACE LOWERING
Surface lowering of a glacier is caused not only by increased ablation, but also by reduced accumulation. The latter factor causes a downstream decrease in ice flux and reduced emergence velocity. Thus, it is essential to estimate the contribution of ablation to surface lowering when considering temporal variations in glacier size. We applied a continuity equation for estimating the contribution of ablation to surface lowering, although only for Area 4, given the availability of flow data for this area (the continuity equation requires historical velocity data). Quincey et al. (2009) measured the flow velocity for the entire Khumbu Glacier, compiling surface velocity maps as several snapshots at 1–2 year intervals for the period 1995–2002. Their maps show the detailed distribution of surface velocity. However, the velocity would have been affected by seasonal variations due to the short interval employed in the study. We did not calculate the continuity equation for Area 3 because of the effect of seasonal variations and the short interval. Hence, we calculated the continuity equation only for Area 4.
The continuity equation expresses conservation of mass for glacier flow with constant density (Paterson, 1994). Consequently, if two of three changes (glacier surface elevation, mass balance, and emergence velocity) are known, the remaining factor can be calculated by the continuity equation. We acquired changes in glacier surface elevation and emergence velocity based on field measurements. However, direct mass balance measurements (e.g., ice stakes) have not been performed. Hence, we adapted the continuity equation for estimating mass balance, which is equivalent to the ablation at Area 4. A change in glacier thickness within a certain part of a glacier can be described by the following continuity equation:
where Q is the ice flux entering and leaving the area of interest; H is ice thickness; t is time; b is the surface mass balance, which is equivalent to the ablation of debris-covered ice;
is the average glacier width (882 m); and x is the longitudinal length of the area of interest (396 m) (Fig. 5). The second term in Equation (2) describes the emergence velocity resulting from glacier flow convergence (Paterson, 1994).
Ice flux at the boundaries of the area of interest is described as follows:
where, W , h, and v are the glacier width, depth, and flow velocity, respectively, at the upper edge of the box for the continuity equation (UECE) and at the lower edge of the box for the continuity equation (LECE) (Fig. 3). In the present study, glacier width (W) was measured from topographical maps (949 m for UECE and 815 m for LECE) and glacier thickness (h) was obtained from the difference between the elevation of the glacier base determined by radio-echo sounding in 1999 (Gades et al., 2000) and the elevation of the glacier surface for each period (345 m for UECE and 328 m for LECE during 1978–1995, and 336 m for UECE and 314 m for LECE during 1995–2004). We applied the surface flow velocities (v) measured by Seko et al. (1998) and those calculated in the present study.
Seko et al. (1998) have estimated the flow velocity data during 1987–1995 using SPOT HRV images (10 m resolution). Unfortunately they have not performed error estimation of those velocities. Considering measurement errors of surface characteristics displacement as at most about 1–2 pixels, those measured displacements could have errors at most 20 m ( = 2.5 m a−1 for flow velocity during 1987–1995). Hence, we considered their estimated values to have small errors and to be reliable.
The flow velocity data measured during 1987–1995 do not strictly correspond to the UECE or LECE values; hence, linear interpolation was performed to obtain the desired values. We modeled the longitudinal distribution of flow velocity using linear regression. Figure 6 shows the longitudinal distribution of flow velocity during 1987–1995. Large residual errors are seen for the upper continuous ice pinnacle zone (distance range in Fig. 6, 0–1000), whereas small residual errors are seen for the separated ice pinnacle zone (distance range in Fig. 6, 1500–2500). The interpolated values are similar to nearby measured values; consequently, we consider the interpolated values to be reliable. In addition, temporal-weighted interpolation was performed to obtain data for the period 1978–1995 from data during 1978–1984 and 1987–1995. We acquired flow data at the upper and lower boundary of the continuity equation. The flow velocity averaged over the depth of the glacier was set as 80% of the mean surface velocity according to Paterson (1994) and Sakai et al. (2006). Other parameters are considered in the following section.
FIGURE 5
Schematic diagram of the volume of glacier considered by the continuity equation. The central cube is the box considered for calculations using the continuity equation. Blue arrows denote glacier flow direction. Transparent boxes on each side of the central cube represent the input and output ice fluxes.
FIGURE 6
Longitudinal distribution of glacier flow velocity during 1987–1995. The data is from Seko et al. (1998). Flow velocities at UECE and LECE are estimated from linear regression of other data.
Results and Discussion
TEMPORAL CHANGES IN ELEVATION
The surface of the investigated parts of the debris-covered Khumbu Glacier showed a significant lowering between 1995 and 2004 (Figs. 7 and 8). The rate of lowering was greater in Areas 2 and 3 than in Areas 1 and 4.
FIGURE 7
DEM differentiation during (a) 1978–1995, and (b) 1995–2004 at 30 m resolution. Thin white line denotes the outline of Khumbu Glacier. The background image is an ASTER image for November 2004.
FIGURE 8
Rate of changes in surface elevation during 1978–1995 and 1995–2004 in this study and those during 1978–1995 by Kadota et al. (2000). Error bars represent the estimated error.
We summed the square of the standard error of elevation changes in each area, as well as the square of the DEM accuracy (Table 1), and presented these values as error bars in Figure 8. Comparison with a previous field study during the period 1978–1995 (Kadota et al., 2000) revealed a significant acceleration of surface lowering in Areas 2 and 3 during the period 1995–2004, especially in Area 3. A minor acceleration was found in Area 1. Areas 2 and 3 are topographically rugged areas that have expanded in size over time (Iwata et al., 2000). Based on remote sensing data, Bolch et al. (2008a) estimated temporal changes in the elevations of glaciers in the Khumbu region with reference to CORONA data (1962), and ASTER data (2002). Their results revealed changes of 0.1–0.6 m a−1 in Area 1 and 0.1–1.3 m a−1 in Areas 2 and 3. These values are slightly smaller than those of the present study (0.4–1.0 m a−1 in Area 1; 0.7–2.3 m a−1 in Area 2; 1.1–2.7 m a−1 in Area 3). However, considering the difference in analysis periods between the previous study (1962–2002) and the present study (1978–1995 and 1995–2004), the changes in elevation appear to be consistent.
A comparison of DEMs reveals the characteristics of temporal changes in elevation at each of four areas in the ablation zone of Khumbu Glacier (Table 2). A relatively small SD of elevation change at Area 1 reflects the slower surface flow speed and smoother topography in this area, indicating greater stability than that in the upper parts of the glacier. The large SDs obtained for Areas 2 and 3 probably reflect the rugged nature of the glacier surface in these areas, where glacier movement resulted in heterogeneous changes in elevation. In fact, Iwata et al. (2000) reported that an area of pronounced relief expanded both upstream and downstream from Area 2 during the period from 1978 to 1995. In this area in particular, the debris cover is not ubiquitous, as there also exist ice cliffs and ponds. Increased surface roughness would have resulted in enhanced heat absorption and thereby increased melting of ice (Sakai et al., 2002). Expansion of the area of pronounced relief from Area 2 to Area 3 would have resulted in accelerated surface lowering in Area 3.
TABLE 2
Statistical summary of the elevation differences among different DEM combinations.
TEMPORAL AND SPATIAL CHANGES IN FLOW VELOCITIES
Table 3 and Figure 9 show spatial and temporal changes in the surface flow velocity of Khumbu Glacier since 1956. Figure 9 shows two trends in the spatial-temporal distribution of surface velocity: a downstream decrease in velocity and a decrease in velocity over time. The former trend is typically observed in glaciers, whereas the latter indicates that the glacier has changed from a steady state to a shrinking state. Other recent studies have also reported decreasing flow velocity at the Khumbu Glacier, based on remote sensing data (Luckman et al., 2007; Quincey et al., 2009).
FIGURE 9
Temporal changes in surface flow speed in Area 4 since 1956. Changes in surface flow speed are modified from Seko et al. (1998), using data from the present study. Vertical axis shows the downstream distance from the base position. Thin arrows correspond to ice gaps G1 and G2, which are located on the upper side of Pinnacle A and C, respectively. Thin arrows for the period 1987–1995 correspond to flow velocities measured by SPOT (Seko et al., 1998). Two dashed lines are the upper and lower edges of the continuity equation (UECE and LECE, respectively), with the upper side representing the input side of ice flux, and the lower side representing the output side. The position of each edge is shown in Figure 3.
TABLE 3
Summary of the surface flow speed (m a−1) for Khumbu Glacier from 1956 to 2004. Flow data are modified from Seko et al. (1998), along with data of the present study. Data for 1956, 1978, and 1984 are from detailed topographical maps (Müller, 1959; Iwata et al., 1980; National Geographic Magazine, 1988) for which the accuracy is unknown. Data for 1987 are from a SPOT HRV image. Data for 1995 and 2004 are based on ground survey data with an accuracy of less than 1 m.
EMERGENCE VELOCITY
We evaluated error propagation from the original errors to emergence velocity and mass balance. We estimated the original errors associated with uncertainty in the glacier width (W), glacier thickness (h), and glacier flow velocity (v) in Equations (2) and (3). The glacier width was measured from a topographic map (yielding values of 949 m for UECE and 815 m for LECE). We estimated the map reading error to be 0.2 mm of the scale factor, following D'Agata and Zanutta (2007). Accordingly, the error in glacier width was calculated to be 20 m. Hence, we estimated the glacier width to be 949 ± 20 m for UECE and 815 ± 20 m for LECE. The elevation of the glacier base (4858 m for UECE and 4855 m for LECE) was measured by radio-echo sounding in 1999 (Gades et al., 2000), with an associated error of 5–20 m. Hence, we estimated the glacier thickness to be 345 ± 20 m for UECE and 328 ± 20 m for LECE during 1978–1995, and 336 ± 20 m for UECE and 314 ± 20 m for LECE during 1995–2004. We evaluated the error in flow velocities during 1978–1995, in terms of the map reading error, following the method proposed by D'Agata and Zanutta (2007). The calculated position errors are 1 m for the 1978 map and 10 m for the 1984 map. These errors propagate to an error in flow velocity for the period 1978–1984 of 41 ± 1.8 m a−1 for both the upper gap (G1) and lower gap (G2) in the pinnacle zone (Table 3). Thus, we hypothesize the following three cases of flow velocity: a downstream decrease in flow velocity (42.8 m a−1 at G1; 39.2 m a−1 at G2), constant flow velocity (41 m a−1 at both G1 and G2), or a downstream increase in flow velocity (39.2 m a−1 at G1; 42.8 m a−1 at G2). It is reasonable to assume that the third case is unrealistic because the rate of glacier flow typically shows a gradual downstream decrease in areas below the equilibrium line altitude (ELA). We estimated the flow velocity at UECE and LECE during 1978–1984 using extrapolation from the flow velocity data at G1 and G2. Subsequently, we performed a temporal weighted interpolation (Table 3).
Finally, we calculated the emergence velocity and mass balance using a continuity equation from the estimated values and the errors outlined above. The emergence velocity and mass balance with errors were 5.72 ± 3.90 m a−1 and ±6.42 ± 3.90 m a−1, respectively. The values between 1995 and 2004, based on ground survey data, were obtained with an accuracy of less than 1 m.
CAUSES OF SURFACE LOWERING
We calculated temporal changes in surface elevation, emergence velocity, and, as the residual term, mass balance, as shown in Table 4. The difference in flow velocity between UECE and LECE has increased in the past decade (3.2 m a−1 during 1995–2004 compared with 2.4 m a−1 during 1978–1995), as shown in Table 3. However, the emergence velocity has decreased due to reduced overall flow speeds arising from a downstream narrowing of the glacier width from 949 to 815 m.
TABLE 4
Temporal changes in elevation, emergence velocity, and surface mass balance in Area 4 during the periods 1978–1995 and 1995–2004.
The residual mass balance, which is equivalent to ablation, has decreased from −6.6 (1978–1995) to −5.66 (1995–2004) m a−1. The estimated emergence velocity have also decreased from 5.9 (1978–1995) to 5.06 (1995–2004) m a−1 (Table 4). Although ablation has shown a decrease, surface lowering has been almost constant due to compensation by reduced emergence velocity. The decrease in ablation at Area 4 seems inconsistent with the present warming climate and with the shrinkage of Himalayan glaciers observed in recent decades (Fujita et al., 1997, 2001). Consequently, the most likely explanation of reduced ablation is enhanced insulation by melt-out debris from the glacier ice, as indicated by the shrinkage of ice pinnacles.
The observed reduction in emergence velocity indicates a decrease in ice flux from the upper catchment. In fact, mountaineers seeking to climb Mt. Everest have reported that the famous icefall dividing the debris-covered ablation area from the accumulation area (e.g., the Western Cwm) has become smoother over the past 50 years. Furthermore, the smoothing of the icefall is apparent in the photographs taken in 1975 and 2004 (Fig. 10). Those support the interpretation of reduced ice flux from upstream parts of the glacier.
FIGURE 10
Photographs of the icefall of Khumbu Glacier in 1975 (left) and 2004 (right). The photograph in 1975 has been provided by the Yomiuri Shinbun Company.
The decrease in ice flux from the icefall may reflect two factors: an increase in melt or a decrease in snowfall. The ELA of Khumbu Glacier is located at the icefall (at around 5700 m a.s.l.; see Benn and Lehmkuhl [2000]), where recent warming is expected to have caused a significant increase in the rate of ice ablation, because the ELA of a glacier is thought to be more sensitive to a change in temperature than are other parts of a glacier reflecting changes in surface albedo (Fujita, 2008a).
It is plausible, therefore, that the present warming climate (Shrestha et al., 1999; Liu and Chen, 2000) has smoothed the icefall via an increase in the rate of ice melt and thus a decrease in ice flux into Area 4. Summer-accumulation-type glaciers are known as sensitive to air temperature (Fujita, 2008a, 2008b). Warming changes the phase of precipitation from solid to liquid. Such a decrease in accumulation can also result in reduced ice flux. Previous studies of Himalayan ice cores have reported a decrease in accumulation on Rongbuk Glacier, on the north side of Mt. Everest (Qin et al., 2002; Kaspari et al., 2008), and on Dasuopu Glacier, on the north side of Mt. Xixabangma (Duan et al., 2006), during the past century.
It is difficult to accurately determine the nature of temporal variations in the dimensions of the accumulation area without field observations. Bolch et al. (2008a) calculated temporal changes in surface elevations in part of the Eastern Khumbu region using two time-separated DEMs (1962 and 2002). The distribution of DEM differentiation shows surface lowering of the accumulation area of some debris-covered glaciers (including Khumbu Glacier). However, the DEMs used in their study include uncertainties in the accumulation area. In contrast, an increase in the size of the accumulation areas of large glaciers, flowing to the south, was reported based on a comparison of maps compiled in the 1950s and 1990s (Salerno et al., 2008). The increase in size was ascribed to the glaciers' favorable orientation in terms of capturing monsoon precipitation. However, their analysis involves uncertainties in terms of interpreting which parts of the basin are glacier, left-over avalanche debris, or disconnected snow patches. Therefore, validation using field measurements is important when analyzing temporal variations in glacier dimensions using remote sensing data. The ongoing discussion on recent changes in the dimensions of the accumulation areas of glaciers requires additional analyses based on field observations or accurate DEMs.
Conclusions
We clarified recent variations in elevation of the debris-covered ablation area of Khumbu Glacier, one of the most well-known debris-covered glaciers in the Himalayas, based on field surveys and analyses of remote sensing data. We assessed the rate of surface lowering at four sites in the debris-covered part of the glacier. Lowering was significant in the middle part of the debris-covered area, consistent with the results of other studies based on remote sensing data, such as Bolch et al. (2008a). Accelerated lowering was found in Areas 2 and 3, in the middle part of the debris-covered area, caused by an increase in the area of topographically rugged glacier surface (Iwata et al., 2000), where ice melt is enhanced at ice cliffs and ponds (Sakai et al., 2000, 2002). To identify the cause of surface lowering at Area 4, in the uppermost part of the debris-covered area, we calculated the contributions to lowering of ablation and emergence velocity. Surface flow velocities have shown a steady decrease since 1956, consistent with the finding of other studies based on remote sensing data (Luckman et al., 2007; Quincey et al., 2009). We found decreases in ablation rate and emergence velocity in the decades before and after 1995. These results suggest that lowering of the glacier surface in the upper ablation area of Khumbu Glacier was caused not only by increasing air temperature, but also by reduced accumulation, although further surveys and analyses of the accumulation area are required to test these conclusions.
Acknowledgments
We wish to thank the staffs of the Department of Hydrology and Meteorology, Ministry of Science and Technology, Nepalese Government, and local sherpas, for their generous assistance during field observations. Valuable comments by K. Nishimura significantly improved the manuscript. A. Sakai helped to improve the manuscript. Field observations and analyses were supported by grants from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan (Nos. 13373066 and 19253001), and by G-4, the 21st Century COE Program of MEXT.
