Introduction
Many of the glaciers in the Tien Shan, Karakoram, Kunlun, Himalaya and Hengduan ranges have extensive mantles of supraglacial debris in their ablation areas (Reference Li and SuLi and Su, 1996; Reference Benn, Wiseman and HandsBenn and others, 2001; Reference Shi, Liu, Wang, Liu and YeShi and others, 2005; Reference Zhang, Liu and DingZhang and others, 2007). Meltwater is an important component of the water cycle in these regions (Reference Yao, Wang, Liu, Pu and LuYao and others, 2004). The wide distribution of supraglacial debris mantles on these glaciers strongly affects the rate of ice melting, which in turn modifies the spatial pattern of mass balance and glacier response to climate change (Reference ØstremØstrem, 1959; Reference McSaveneyMcSaveney, 1975; Reference Nakawo and YoungNakawo and Young, 1981; Reference Mattson, Gardner and YoungMattson and others, 1993; Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Nicholson and BennNicholson and Benn, 2006). An empirical relationship between supraglacial debris thickness and ice melt rate has been established based on experimental observations of glaciers (Reference ØstremØstrem, 1959; Reference Nakawo and YoungNakawo and Young, 1981; Reference Mattson, Gardner and YoungMattson and others, 1993; Reference Adhikary, Seko, Nakawo, Ageta and MiyazakiAdhikary and others, 1997; Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000) and laboratory experiments (Reference Reznichenko, Davies, Shulmeister and McSaveneyReznichenko and others, 2010), which reveals that debris thickness principally controls the sub-debris ablation rate. Therefore, estimation of the extent and thickness of debris cover and its effect on ice ablation is crucial in determining glacier runoff and evaluating water resources in affected regions.
Because of practical difficulties in surveying on debris-covered glaciers, detailed in situ measurements of debris thickness have only been performed on a few debris-covered glaciers (e.g. Reference Nakawo, Iwata, Watanabe and YoshidaNakawo and others, 1986; Reference Nicholson and BennNicholson and Benn, 2006; Reference MihalceaMihalcea and others, 2008a,Reference Mihalceab) and estimating the spatial variation of debris thickness with sufficient accuracy has thus far proven difficult (Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997; Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007; Reference Racoviteanu, Paul, Raup, Khalsa and ArmstrongRacoviteanu and others, 2009). Currently, satellite data have the potential to provide information over wide spatial and temporal scales for monitoring glaciers across large and remote glacierized regions (Reference Raup, Racoviteanu, Khalsa, Helm, Armstrong and ArnaudRaup and others, 2007; Reference Racoviteanu, Paul, Raup, Khalsa and ArmstrongRacoviteanu and others, 2009). Several approaches to debris-cover mapping have been applied to different debris-covered glaciers (e.g. Reference Bishop, Kargel, Kieffer, MacKinnon, Raup and ShroderBishop and others, 2001; Reference Taschner and RanziTaschner and Ranzi, 2002; Reference Paul, Huggel and KääbPaul and others, 2004; Reference Racoviteanu, Williams and BarryRacoviteanu and others, 2008). However, there is no single best method for debris-cover mapping that can be applied to large regions without some manual corrections of the resulting outlines; these various methods have not yet been compared and a superior method has thus not yet emerged (Reference Racoviteanu, Paul, Raup, Khalsa and ArmstrongRacoviteanu and others, 2009). Some investigations of debris-covered glaciers using satellite data have focused on extracting surface temperature information for numerical models of ice melt beneath debris cover (Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997; Reference Nakawo and RanaNakawo and Rana, 1999) or deriving supraglacial debris-cover patterns (Reference MihalceaMihalcea and others, 2008b). To assess the effect of debris cover on ice ablation, several numerical models have been proposed, generally based on energy balance for the debris layer and driven by meteorological variables and the physical properties of debris covers (e.g. Reference McSaveneyMcSaveney, 1975; Reference Nakawo and YoungNakawo and Young, 1981, Reference Nakawo and Young1982; Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Han, Ding and LiuHan and others, 2006; Reference Nicholson and BennNicholson and Benn, 2006; Reference Reid and BrockReid and Brock, 2010). The main challenge in applying these models is the requirement for high-quality input parameters related to the extent, thickness and thermal properties of debris covers. Field determinations of such parameters are especially problematic on debris-covered glaciers. For this reason, the parameter of ‘thermal resistance’ was proposed, defined as the debris thickness divided by the thermal conductivity of the debris layer (Reference Nakawo and YoungNakawo and Young, 1981, Reference Nakawo and Young1982). Some studies have been conducted to derive the distribution of thermal resistance from satellite remote-sensing data (Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997; Reference Nakawo and RanaNakawo and Rana, 1999; Reference Suzuki, Fujita and AgetaSuzuki and others, 2007). However, the representativeness of thermal resistance as debris thickness has not been confirmed and the performance of a surface energy-balance model that accounts for the thermal resistance of the debris layer has only been validated with hydrographs.
This study had the following goals. First, we investigated the spatial distribution of debris thickness throughout the ablation area of Hailuogou glacier, southeastern Tibetan Plateau, based on detailed in situ survey data and recent satellite images. High-resolution in situ survey measurements of debris thickness were compared with the spatial pattern of thermal resistance derived from satellite images to assess the validity of applying thermal resistance as a proxy for debris-thickness distribution. Second, to understand the role of debris thickness in spatial variability of melt rate, the ice melt rate of the entire ablation area was estimated by a surface energy-balance model that took thermal resistance of the debris layer into account and was forced using daily mean meteorological variables. Model performance was evaluated by comparing calculated melt rates with sub-debris melt rates measured by ablation stakes. The findings highlight the role of spatial variability in debris thickness in the spatial characteristics of ice melt and its correlated processes on a debris-covered maritime glacier. Such work is a necessary first step toward understanding mass-balance variation of the debris-covered maritime glacier and its response to climate changes in the Tibetan Plateau.
Study Area
Hailuogou glacier is the longest maritime glacier located on the eastern side of Mount Gongga (29°35.6′ N, 101°56.7′ E; Fig. 1). The glacier has an area of ∼25.7 km2, a length of 13.1 km and an altitude range of 2901–7556 m a.s.l. (Reference Shi, Liu, Wang, Liu and YeShi and others, 2005), generally flowing eastward. Below an icefall (3650–4980 m a.s.l.; Fig. 1), debris covers most of the glacier tongue due to high rates of debris supply from mixed ice/snow/rock avalanching from the icefall and surrounding rock walls through frost weathering processes and structural rockfalls (Reference Li and SuLi and Su, 1996). The surface debris consists primarily of granite, metasandstone and slate dust, silts, sands, gravels, rocks and boulders, with primary mineral components including quartz, potassium feldspars (in the form of orthoclase and microcline), feldspars (albite, calcium-rich albite), micas (muscovite, biotite, phlogopite) and carbonates (calcite) (Reference Li and SuLi and Su, 1996). The dominant climatic influences on the glacier are the southeast monsoon in summer and westerly circulation in winter (Reference Li and SuLi and Su, 1996; Reference Su, Song and CaoSu and others, 1996). The glacier receives abundant average annual precipitation (average value of ∼1.9 m w.e. for 1988–2004 at the terminus) of which summer precipitation (May–October) accounts for ∼80% of the annual total. The average annual temperature at the terminus was 4.1°C over the period 1988–2004 (Reference Zhang, Fujita, Liu, Liu and WangZhang and others, 2010).
Fig. 1. Map of the ablation area of Hailuogou glacier on the eastern side of Mount Gongga, southeastern Tibetan Plateau. The colour image is a false-colour Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) image acquired on 18 January 2009. The location of the automatic weather station (AWS) is marked. Dashed lines (P1–P5) are transverse and central line profiles along which debris thickness was measured. GAEORS: Gongga Alpine Ecosystem Observation and Research Station.
Reference Zhang, Fujita, Liu, Liu and WangZhang and others (2010) assessed the nature of multidecadal changes in surface elevation and ice velocity below the icefall, based on digital elevation models created for different periods and in situ differential GPS survey data. The authors determined that the glacier below the icefall thinned by 1.1 ± 0.4 m a−1 during 1966–2009 and the average summer ice velocity in 2008 had undergone a significant reduction compared with that of the 1980s. The terminus receded at a rate of 12.7 m a−1 during 1966–89 (Reference Su, Liu, Wang and ShiSu and others, 1992) and at a rate of 27.4 m a−1 during 1998– 2008 (Reference Zhang, Fujita, Liu, Liu and WangZhang and others, 2010).
Data
In situ survey data
Extensive in situ surveys of debris thickness were conducted during the summer of 2009 across the altitudinal range 2900–3600 m a.s.l. A total of 300 points were surveyed over the entire ablation area. At each measurement site, we excavated the debris cover until reaching glacier ice (Fig. 2). We then laid a straight edge across the top of the hole and measured the vertical distance from the ice surface to the top of the debris cover (at the straight edge). If the supraglacial debris was <20 cm thick, at least three holes were measured at each site. Using these data, transverse and longitudinal profiles were constructed across representative areas of the glacier (Fig. 1), each based on eight or more samples obtained by excavation.
Fig. 2. Distribution of debris thickness at different elevations in the ablation area of Hailuogou glacier. Inset shows the method of in situ surveys of debris thickness in the ablation area during the summer of 2009.
Ablation data were collected during the summer of 2008 using 27 ablation stakes drilled into the glacier at different elevations (Fig. 1). The stakes were arranged in one central longitudinal profile (mainly east–west) and six transverse profiles (north–south). Readings were taken at 7 day intervals.
An automatic weather station (AWS) operated in the upper region of the ablation area at ∼3550 m a.s.l. (Fig. 1). At this site, the surface is relatively flat and homogeneous. The AWS was installed in April 2008, and the data were analyzed for the period 13 May 2008 to 9 August 2009. Sensor specifications are listed in Table 1. Air temperature, relative humidity, wind speed and wind direction were measured at 2.0 and 4.0 m above the glacier surface. A radiation sensor measured the incoming and reflected solar radiation at 1.5 m above the glacier surface. A Campbell Scientific CR10× data logger stored 10 min averages of all the above variables. This study also used precipitation and temperature data from a precipitation gauge T-200 located at the glacier margin (3500 m a.s.l.; Fig. 1) for the period 1 October 2008 to 9 August 2009 and from a sub-alpine observatory (Gongga Alpine Ecosystem Observation and Research Station, GAEORS) for the period 13 May to 30 September 2008, 1.5 km from the terminus (3000 m a.s.l.; Fig. 1). All variables were compiled to daily values.
Table 1. Specifications of instruments at the automatic weather station operated in the upper region of the ablation area of Hailuogou glacier
ASTER data
Two orthorectified imageries of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) on NASA’s Terra satellite with a spatial resolution of 90 m were used for this analysis, obtained on 17 December 2008 and 18 January 2009 (1158 h local time). The imageries were generated and distributed as Level 3A01 products by the ASTER Ground Data System (ASTER GDS) at the Earth Remote Sensing Data Analysis Center (ERSDAC) in Japan (Reference Fujisada, Bailey, Kelly, Hara and AbramsFujisada and others, 2005; www.gds.aster.ersdac.or.jp/gds_www2002/exhibition_e/a_products_e/a_product2_e.html). The location of the images was affine-transformed by referring to a topographical map, and the root-mean-square error (rmse) of the affine transformation was within 15m (Reference Zhang, Fujita, Liu, Liu and WangZhang and others, 2010).
Methods
Energy-balance model
The energy balance at the debris surface is expressed as
where R n, H, LE, P and G are net radiation flux, net sensible heat flux, net latent heat flux, heat flux by rain on the surface and conductive heat flux into the debris layer, respectively. All terms were taken to be positive towards the debris surface (W m−2).
Net radiation
The net radiation flux, which consists of short- and longwave radiations, is given by
where R S is downward shortwave radiation, R L is downward longwave radiation, α is albedo, ε is the emissivity of the debris surface (taken to be 1), ρ is the Stefan–Boltzmann constant (5.67 × 10−8 W m−2 K−4) and Ts is the surface temperature of the debris (°C). Here, upward longwave radiation was calculated from the Stefan–Boltzmann law and downward longwave radiation was calculated from air temperature, relative humidity and the ratio of downward shortwave radiation to that at the top of the atmosphere using an empirical scheme (Reference KondoKondo, 1994).
Turbulent heat fluxes
Turbulent sensible and latent heat fluxes are calculated using the bulk method, expressed as follows:
where c a is the specific heat capacity of air (1006 J kg−1 K−1), ρ a is the density of air (kg m−3), C is the bulk coefficient for sensible and latent heat (0.002), U is wind speed (m s−1), Ta is air temperature (°C), rh is relative humidity, l e is the latent heat of the evaporation of water (2.50 × 106 J kg−1) and q is saturated specific humidity (kg kg−1). This latter value was calculated as a function of air temperature using the empirical equations presented by Reference KondoKondo (1994):
where p a and e are air and vapour pressure (Pa) respectively.
Heat flux supplied by rain
The heat flux supplied by rain is often ignored in debris-covered glacier studies because it tends to be extremely small in comparison with other heat fluxes. Therefore, relatively few studies have considered the effect of rain on ice melt rates (e.g. Reference McSaveneyMcSaveney, 1975; Reference Sakai, Fujita and KubotaSakai and others, 2004; Reference Hock and HolmgrenHock and Holmgren, 2005; Reference Reid and BrockReid and Brock, 2010; Reference Reznichenko, Davies, Shulmeister and McSaveneyReznichenko and others, 2010). Here the energy supplied by rain, P, is calculated by an equation suggested by Reference Hay and FitzharrisHay and Fitzharris (1988) and Reference Hock and HolmgrenHock and Holmgren (2005):
where ρ w is the density of water (999.7 kg m−3), c w is the specific heat capacity of water (4181.3 J kg −1 K−1), w is the rainfall rate in m s−1 and Tr is the rainfall temperature (°C; assumed to be equal to air temperature).
Calculating sub-debris ice melt rate
The only heat flux considered to reach the glacier ice is the conductive heat flux with the simplifying assumption of a linear temperature profile within the debris layer and the constant heat flux stored in the debris layer from day to day (e.g. Reference KrausKraus, 1975; Reference Nakawo and YoungNakawo and Young, 1981; Reference Nicholson and BennNicholson and Benn, 2006). Therefore, the heat flux used for ice ablation, Q m, can be calculated as
where R is thermal resistance of the debris layer (m2 K W−1) and Ti is the temperature at the ice-debris surface (°C). Previous investigations on Hailuogou glacier showed that the ice temperature is approximately 0°C at the ice-debris surface (Reference Li and SuLi and Su, 1996; Reference Su, Song and CaoSu and others, 1996).
Following from the above, the sub-debris melt rate, M (m w.e.d−1), was calculated as
where ρ i is the density of ice (900 kg m−3) and L f is the latent heat of fusion (3.34 × 10−5 J kg−1).
As shown in Equations (2–7), all energy-balance components, except the shortwave radiation term, are determined explicitly from the debris surface temperature, Ts. The debris surface temperature can be calculated numerically by iteration, a technique often used in similar studies for different debris-covered glaciers (Reference Nicholson and BennNicholson and Benn, 2006; Reference Reid and BrockReid and Brock, 2010). The resulting daily mean surface temperature is then used to calculate the rate of ice melt from Equations (7) and (8) for prescribed thermal resistance of the debris layer.
To solve the system of Equations (1–8), our model requires albedo, thermal resistance and daily mean meteorological variables (air temperature, relative humidity, wind speed and solar radiation). The surface albedo and thermal resistance of the debris layer were estimated using ASTER data, while the latter variables were 24 hour mean values measured in the summer of 2008.
Calculating bare-ice melt rate
To quantify the effect of the spatial distribution of debris thickness, the bare-ice melt rate, M’ (mw.e.d−1), can be calculated as
where all heat fluxes are calculated using the same equations described above and the downward conductive flux is considered negligible. The albedo of bare ice was derived from ASTER data obtained on 18 January 2009, set to 0.3.
Thermal resistance of a debris layer
Thermal resistance of a debris layer, which is defined as debris thickness divided by the thermal conductivity of the debris layer (Reference Nakawo and YoungNakawo and Young, 1981, Reference Nakawo and Young1982), is expressed as
where H and λ are debris thickness and the thermal conductivity of the debris layer, respectively.
Note that field determinations of the thickness and thermal conductivity of the debris layer are especially time-consuming and not practical on glaciers. Therefore, thermal resistance of debris layers on Hailuogou glacier was calculated from the surface temperature derived from the thermal infrared (TIR) bands of ASTER and the net radiation estimated from the visible/near-infrared (VNIR) bands of ASTER and US National Centers for Environmental Prediction (NCEP)/US National Center for Atmospheric Research (NCAR) reanalysis data that corresponded to the nearest time and location of ASTER acquisition. A schematic diagram of the method is shown in Figure 3b. As can be seen in Figure 3b, the average brightness temperature was calculated from five TIR bands of ASTER and used as the surface temperature on the glacier. Spectral reflectance at the top of the atmosphere in VNIR bands of ASTER was averaged to estimate broadband albedo. The ASTER-derived surface temperature and albedo associated with surface downward radiation fluxes from NCEP/NCAR reanalysis data (Reference KalnayKalnay and others, 1996) that corresponded to the nearest time and location of ASTER acquisition were used to calculate thermal resistances of debris layers on the glacier (Fig. 3). The pixel size of the ASTER image is 90 m and thermal resistance was calculated using the same resolution.
Fig. 3. (a) ASTER VNIR band imagery on 18 January 2009. (b) Schematic diagram of the method of calculating thermal resistance.
This approach is based on the fact that the net radiation is usually the dominant heat source on Himalayan debris-covered glaciers and the contribution of the turbulent heat fluxes to the total energy balance can be neglected (Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Takeuchi, Kayastha and NakawoTakeuchi and others, 2000; Reference Suzuki, Fujita and AgetaSuzuki and others, 2007). Also, we found the same characteristic of the energy balance at the debris-covered surface on Hailuogou glacier. The uncertainties of this approach, which neglects the turbulent heat fluxes in the energy-balance calculation and uses multitemporal ASTER data, were evaluated by Reference Suzuki, Fujita and AgetaSuzuki and others (2007). They found that these uncertainties are unlikely to affect the spatial pattern of thermal resistance. Although shading of solar insolation by surrounding mountains may not affect thermal resistance of the debris layer, the presence of water in the debris layer affects thermal resistances (Reference Suzuki, Fujita and AgetaSuzuki and others, 2007). Therefore, two ASTER imageries were used on 17 December 2008 and 18 January 2009 (1158 h local time). Note that precipitation from May to October on Hailuogou glacier accounts for ∼80% of the annual total (Reference Zhang, Fujita, Liu, Liu and WangZhang and others, 2010), i.e. only a small amount of precipitation falls in winter (November–April). According to the meteorological data observed at the AWS and precipitation gauge (Table 2), there was clear sky without precipitation at the time of ASTER acquisition. Therefore, the debris surface could be considered to be dry in the calculation of thermal resistance.
Table 2. Meteorological conditions at 1200 h on 17 December 2008 and 18 January 2009. T, rh, U and P are the air temperature, relative humidity, wind speed and precipitation, respectively
Previous studies have found that the thermal resistance of the debris layer at a specific site can usually be regarded as constant in a melting season on debris-covered glaciers (Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997; Reference Nakawo and RanaNakawo and Rana, 1999; Reference Suzuki, Fujita and AgetaSuzuki and others, 2007). ASTER-derived thermal resistances are available for 2008 and for 2009 for comparison at the same pixels on Hailuogou glacier. We found that although ASTER-derived thermal resistances are instantaneous, thermal resistances derived from the independent data are not significantly different. The overall correlation coefficient between thermal resistances in 2008 and those in 2009 is 0.88 (significance level p < 0.001). This finding is in agreement with Reference Suzuki, Fujita and AgetaSuzuki and others (2007), who found that there is little difference in thermal resistances derived from the independent data. Thermal resistances of debris layers in the ablation area of Hailuogou glacier were therefore considered as constant in the melting season of 2008.
Results and Discussion
Spatial distribution of debris thickness
Figure 2 shows the distribution of debris thickness at different elevations in the ablation area of Hailuogou glacier. The general trend is to increasing debris thickness with decreasing elevation (Fig. 2), although in some areas the spatial distribution is inhomogeneous because of crevasses, supraglacial ponds and ice cliffs. In situ survey data indicate that debris thickness increases from just several millimetres of patchy cover below the icefall to >1 m at the terminus; in some areas where piles of large rocks exist, the thickness reaches several metres.
The debris cover of Hailuogou glacier is clearly recognized in the ASTER VNIR band imagery (Fig. 3), which appears in grey in the imagery and is distinct from surrounding red-coloured vegetation. Figure 4 shows spatial variations in debris thickness and ASTER-derived thermal resistance along transverse profiles (Fig. 1) that are representative of the main sectors of the glacier (terminus, and central and upper parts of the ablation area) and the longitudinal profile (Fig. 1) that records continuous central flowline variation. The variation in debris thickness is pronounced in profiles P1 and P2, whereas it is subdued in profiles P3 and P4. In profile P1, debris is thicker at the southern and northern lateral margins of the glacier than in the central part. Across profiles P2–P4, debris thickness shows a gradual increase with increasing distance from the southern lateral margin of the glacier (Fig. 4).
Fig. 4. Transverse and longitudinal profiles of debris thickness (points) and ASTER-derived thermal resistance (lines): (a) P1 (3000 m a.s.l.), (b) P2 (3150 m a.s.l.), (c) P3 (3380 m a.s.l.), (d) P4 (3550 m a.s.l.) and (e) P5 (see Fig. 1 for profile locations). S and N denote the southern and northern parts of the ablation area, respectively. The dashed line represents the northern boundary of the glacier.
Detailed field surveys of debris thickness on glaciers are particularly time-consuming and challenging. Therefore, we attempted to use the thermal resistance of the debris layer as a proxy for debris thickness on Hailuogou glacier. In Figure 5a the relation between debris thickness and ASTER-derived thermal resistance is shown. Debris thickness correlated well with thermal resistance, yielding a correlation coefficient of r = 0.61 (p < 0.001), although some scatter occurs in the data (Fig. 5a). The transversal profiles (P1–P4) showed a general trend in which ASTER-derived thermal resistance increased toward the northern lateral margin, which corresponded well to the transverse distribution of known debris thickness on the glacier. Across profiles P1 and P2, the variation of thermal resistance is pronounced due to the inhomogeneous distribution of debris thickness, which varies from a few centimetres to >1 m (Fig. 4a and b). In contrast, debris thickness variation is small in profiles P3 and P4, which corresponded well to the variation of thermal resistance (Fig. 4c and d). Along a central flowline (P5; Fig. 1), thermal resistance showed a gradual increase down-glacier, consistent with the surveyed debris thickness; these variations reflected the down-glacier transition from patchy thin debris to thick continuous debris cover (Fig. 4e). Overall, ASTER-derived thermal resistance reflected the known broad-scale variations in debris thickness.
Fig. 5. (a) Relationship between surveyed debris thickness and ASTER-derived thermal resistance, (b) comparison of modelled and observed ice melt rates and (c) modelled and observed melt rates plotted against debris thickness for Hailuogou glacier. If an ASTER gridcell contained more than one survey point, the average debris thickness is used in (a). Error bars in (a) denote standard errors.
The spatial distribution of ASTER-derived thermal resistance is illustrated in Figure 6a. Relatively high values of thermal resistance were found at the northern lateral margin of the glacier (∼3000–3300 m a.s.l.) where the debris cover is thick and continuous; lower values occurred mainly below the icefall, where the debris cover is thin and discontinuous (Fig. 6a). A low value of thermal resistance was also found at the southern lateral margin (∼3000– 3200ma.s.l.), which is attributed mainly to the complex pattern of surface types in this area (e.g. ice cliffs, supraglacial ponds, bare ice surface and debris-covered surface). Comparison of ASTER-derived thermal resistances along the transverse profiles (Fig. 1) on and outside the glacier shows that thermal resistances of debris layers on the glacier are less than those outside the glacier (Table 3). The difference in thermal resistance on (4.2 × 10−2 m2 K W−1) and outside (5.4 × 10−2 m2 K W−1) the glacier is small in profiles P1 and P2. In contrast, the difference on (1.6 × 10−2 m2 K W−1) and outside (5.7 × 10−2 m2 K W−1) the glacier is pronounced in profiles P3 and P4.
Table 3. ASTER-derived thermal resistances along the transverse profiles on and outside the glacier. Unit is 10−2 m2 K W–1
Fig. 6. Spatial distributions of (a) ASTER-derived thermal resistance (TR) value, (b) melt rate calculated by the presented model and (c) melt ratio defined as the sub-debris melt rate divided by the bare-ice melt rate at the same elevation in the ablation area of Hailuogou glacier.
As noted above, we first compared the spatial distribution of debris thickness and the spatial pattern of ASTER-derived thermal resistance. Overall, the thermal resistance performed well in representing the general pattern of spatial variations in debris thickness on the glacier.
Role of debris thickness in ice-melting pattern
The average value of debris thickness is 0.16 m in the ablation area, based on in situ debris thickness survey data during the summer of 2009. The debris-covered area accounts for 41% of the ablation area on Hailuogou glacier. To understand the role of debris thickness and its spatial continuity in modifying the spatial characteristics of ice melt rates, we present a surface energy-balance model for calculating sub-debris melt from daily mean meteorological variables (Equations (1–8)) in which thermal resistance of the debris layer is taken into account. Model performance was evaluated by comparing modelled melt rates with observations from ablation stakes. Observed and modelled melt rates are compared in Figure 5b. The overall correlation coefficient between observed and modelled melt rate at ablation stakes was 0.82 (p < 0.001), although some points were plotted away from the slope = 1 line. A comparison of modelled and observed melt rates at different thicknesses of debris layer (Fig. 5c) revealed that simulations at different debris thicknesses were in good agreement with observed melt rates. This result confirms that a model that accounts for all components of the energy fluxes and thermal resistance of the debris layer can produce reliable estimates of ice melt rates beneath various debris thicknesses on an annual timescale.
In applying the model to the entire ablation area, the spatial distribution of the melt rate was derived from daily mean meteorological variables and ASTER-derived thermal resistance for the period 13 May to 5 October 2008 (Fig. 6b). The spatial resolution of the calculated melt rate distribution was defined by a 90 m ASTER pixel size. The calculated melt rate in the ablation area ranged from a minimum of 0.023 m d−1 to a maximum of 0.058 m d−1. Low melt rates were seen mainly at the northern lateral margin (3000–3300 m a.s.l.), where the debris cover is thick and continuous, thereby insulating and preserving the underlying ice (Fig. 6b). High melt rates were found in two areas (Fig. 6b): (1) below the icefall (∼3450–3600 m a.s.l.), where the debris cover is thin and patchy; and (2) at the southern lateral margin (3000–3200 m a.s.l.), an area with a wide distribution of coexisting large ice cliffs, supraglacial ponds, bare ice and debris cover. Ice cliffs are generally covered with only a thin layer of debris or dust, meaning that the melt rate is higher than on nearby debris-covered surfaces (Reference Iwata, Watanabe and FushimiIwata and others, 1980; Reference Sakai, Takeuchi, Fujita and NakawoSakai and others, 2000). In addition, the energy absorption at the supraglacial ponds is several times larger than at the surrounding debris-covered surface, resulting in accelerating ice melt (Reference Sakai, Takeuchi, Fujita and NakawoSakai and others, 2000). Hence, these parts of the glacier are likely to make disproportionately large contributions to ablation; meanwhile, the debris layer becomes increasingly unstable due to the heterogeneous ablation in these parts of the glacier, releasing debris via occasional sluicing of debris down steep slopes. This process redistributes the debris upon the glacier surface (Reference Iwata, Aoki, Kadota, Seko and YaaguchiIwata and others, 2000; Reference Benn, Wiseman and HandsBenn and others, 2001) and then modifies the spatial variability in the ice melt rate.
The altitude distribution of debris thickness (Fig. 2) also complicated the ablation gradient on the ablation area, making it completely different to that of debris-free ice (Fig. 7). The mean ablation gradient for debris-free ice was 0.25 m (100 m)−1 calculated using Equation (9), assuming an entirely debris-free surface on the ablation area, whereas the ablation gradient for the real surface was reversed at 0.9 m (100 m)−1 below 3200 m a.s.l. and +0.6 m (100 m)−1 above 3200 m a.s.l. (Fig. 7). The melt ratio was defined as the sub-debris melt rate, M, divided by the bare-ice melt rate at the same elevation, M’, calculated by Equation (9). A melt ratio of 1.0 indicates that the ice melt rate beneath the debris layer equals that of bare ice, whereas a ratio >1.0 indicates that debris enhances the ice melt rate and a ratio <1.0 indicates that debris inhibits the ice melt. The spatial distribution of the melt ratio clearly shows the crucial effect of debris thickness and its spatial continuity on ice melting in the ablation area (Fig. 6c). About 67% of the ablation area has undergone accelerated melting, whereas about 19% of the ablation area has experienced inhibited melting. On only 14% of the ablation area does the sub-debris melt rate equal the bare-ice melt rate (Fig. 6c). The heterogeneous distribution of debris thickness results in marked spatial variability in the ablation regime on Hailuogou glacier.
Fig. 7. Variation in the ablation gradient at debris-covered and debris-free ice for Hailuogou glacier.
Discussion
In the case of debris-covered glaciers, the main challenge for glacial modelling is to acquire information on the extent, thickness and thermal properties of debris covers, especially at a regional scale. Reference Reid and BrockReid and Brock (2010) (Fig. 6) reported that variation in debris thermal conductivity had a significant effect on the performance of their model and pointed out that although the average thermal conductivity for a debris layer can be estimated by observing data, it is difficult to quantify how it may vary with debris depth. Other studies have made similar findings (e.g. Reference Han, Ding and LiuHan and others, 2006; Reference Nicholson and BennNicholson and Benn, 2006). In this study, the thermal resistance derived from ASTER data correlated well with ground-surveyed debris thickness (Fig. 5a), and a comparison of the distribution of supraglacial debris thickness and the spatial pattern of ASTER-derived thermal resistance validated the use of satellite-derived thermal resistance as a proxy for debris thickness over large areas (e.g. 90 m × 90 m). In addition, comparison of ASTER-derived thermal resistance in 2008 with that in 2009 at the same pixels indicates clearly that thermal resistances derived from the independent data are not significantly different. Therefore, the model described by Equations (1–8) calculates the ice melt rate beneath the debris directly from meteorological data and satellite-derived thermal resistance, without the need for the extent, thickness and thermal properties of debris covers, which possibly provides important insight into studying debris-covered glaciers.
The model described above is run with meteorological inputs from AWS sites. We ran it with measurements from just one site and, with the exception of air temperature, have not extrapolated to other areas. There may be an unquantifiable error caused by the above given that variations in altitude, aspect and shading result in very different meteorological data conditions. In addition, the water content of the debris layer may influence the model performance, which is difficult to determine for the entire debris-covered area due to the inhomogeneous structure of the debris layer (with dry rubble and boulders on top containing considerable air spaces, and smaller rubble and fine particles, which may often be saturated with water near the surface). The water content of the debris layer is unlikely to be known in any practical application of the model presented in previous studies (e.g. Reference Nakawo and YoungNakawo and Young, 1981; Reference Nicholson and BennNicholson and Benn, 2006; Reference Brock, Mihalcea, Kirkbride, Diolaiuti, Cutler and SmiragliaBrock and others, 2010; Reference Reid and BrockReid and Brock, 2010). Reference Nicholson and BennNicholson and Benn (2006) tackled the unknown nature of the water content of the debris by considering the two limiting cases: the debris cover being either completely dry or fully saturated. They found that modelled ice melt rates for dry debris were less than those for wet debris. This assumption is undoubtedly not the case in the real world. Our study does not take into account the effect of the water content of the debris layer since considering the effect of the water content must allow the physical properties of the debris layer to vary with depth, which would require detailed data on internal debris temperature and water content. This effect may explain why modelled melt rates were less sensitive to debris thickness than observations (Fig. 5c). This phenomenon should be a focus of future studies.
Note that ice ablation is controlled by the amount of energy supplied to the ice (Reference PatersonPaterson, 1994), which is controlled in turn by the interplay between the surface energy balance (which determines energy received at the surface) and heat transfer processes through the debris layer (which determine the amount of energy reaching the ice interface) (Reference PatersonPaterson, 1994; Reference Benn, Wiseman and HandsBenn and others, 2001). Figure 8a shows trends in the surface energy-balance components with increasing debris thickness, as calculated at five ablation stakes around 3150 m a.s.l. using the model described above. Overall, the net radiation is the main heat source at the debris surface and it is stable with increasing debris thickness; secondary heat sources are turbulent sensible and latent heat fluxes, which show a gradual decrease with increasing debris thickness (Fig. 8a). The conductive heat flux into the debris is a significant heat sink, decreasing with increasing debris thickness (Fig. 8a). Although Hailuogou glacier receives abundant precipitation in the summer, the heat flux supported by rain has little effect on the surface energy balance at the debris surface (Fig. 8a), accounting for only 2% of the total energy components. However, such results should be treated with caution because precipitation data were not available for the glacier during the summer of 2008 and the data used were from a sub-alpine observatory (GAEORS) 1.5 km from the terminus (3000 m a.s.l.; Fig. 1). For this reason, we calculated the ice melt rate of the entire ablation area over the period 1 May to 9 August 2009, when precipitation data on the glacier were available. As shown in Figure 9, there is no significant difference between the ice melt rates considering the effect of rain and those neglecting the effect of rain. The characteristic of the debris surface energy balance on Hailuogou glacier (Fig. 8) is in agreement with the findings of previous studies in the Himalaya (Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Takeuchi, Kayastha and NakawoTakeuchi and others, 2000; Reference Suzuki, Fujita and AgetaSuzuki and others, 2007), the Karakoram mountains (Reference Mattson and GardnerMattson and Gardner, 1989), the Italian Alps (Reference Nicholson and BennNicholson and Benn, 2006; Reference Brock, Mihalcea, Kirkbride, Diolaiuti, Cutler and SmiragliaBrock and others, 2010) and Svalbard (Reference Nicholson and BennNicholson and Benn, 2006). These studies suggested that the debris surface energy balance is dominated by the pattern of net radiation which heats the debris, and the contribution of the turbulent heat fluxes to the total energy balance is normally small; the heat flux from precipitation is normally ignored as it is negligibly small compared with the other fluxes. Reference Reid and BrockReid and Brock (2010) made a statistical analysis of the effect of precipitation at a debris surface and suggested that the heat flux from precipitation can be omitted without large losses in model performance, which is in agreement with our results (Fig. 9).
Fig. 8. (a) Relationship between surface energy-balance components and debris thickness at a given elevation and (b) relationship between surface energy-balance components and elevation for a debris layer with a thickness of 0.02 m for Hailuogou glacier. R n, H, LE, P and G are net radiation flux, net sensible heat flux, net latent heat flux, heat flux by rain on the surface and conductive heat flux into the debris layer, respectively.
Fig. 9. Comparison of the ice melt rate calculated by considering and neglecting the effect of rain for Hailuogou glacier.
At the same altitude, an increase in debris cover from 0.02 m to 0.13 m is accompanied by a 36% decrease in the energy available for melting (Fig. 8a), showing that a small change in debris thickness results in a marked change in the energy available for ice melting. It is also important to consider the effect of decreasing energy input to the debris layer with increasing elevation, as indicated by the different energy components obtained for areas with a constant debris cover of 0.02 m but located at different elevations (Fig. 8b). This effect may explain why we were unable to determine the effective thickness (where the sub-debris melt rate is maximized) and the critical thickness (where the sub-debris melt rate is equal to the bare-ice melt rate, and with increase of which the debris has an insulation effect on the underlying ice) from observed ablation data for the summer of 2008.
Conclusions
Supraglacial debris thickness is a key element in estimating the energy and mass balances for Hailuogou glacier. In situ surveys reveal that debris thickness increases from several millimetres of patchy cover below the icefall to >1 m at the terminus. An analysis of transverse and longitudinal profiles demonstrates a strong relationship between debris thickness and ASTER-derived thermal resistance. The spatial pattern of thermal resistance reflects broad-scale variations in debris extent and thickness, thereby supporting the use of thermal resistance as a proxy for debris thickness over large areas.
A surface energy-balance model has been presented, involving all components of the energy fluxes at the debris surface, which is able to accurately produce ice melt rates beneath various debris thicknesses on an annual timescale. By accounting for thermal resistance of the debris layer, this approach calculates the ice melt rate beneath the debris directly from meteorological data, without the need for detailed field surveys of the extent, thickness and thermal properties of debris covers. The main advantage of this approach is the requirement for fewer input variables and parameters. The data used, namely meteorological data and ASTER data, are likely to be available. Even if there are no meteorological observations, these data can come from reanalysis. Hence, this approach can be applied to the majority of debris-covered glaciers.
The ice-melt rates beneath debris covers on the entire ablation area of Hailuogou glacier were calculated by the model from daily mean meteorological variables and ASTER-derived thermal resistances. The nature of spatial variability in the ice melt rate reveals that debris thickness and its spatial continuity exert a strong influence on the spatial distribution of glacier melting: about 67% of the ablation area has undergone accelerated melting, whereas about 19% of the ablation area has experienced inhibited melting, and the sub-debris melt rate equals the bare-ice melt rate in only 14% of the ablation area.
Acknowledgements
We thank the Gongga Alpine Ecosystem Observation and Research Station of the Chinese Ecological Research Network for providing the meteorological data. This research was supported by the National Basic Research Program of China (2007CB411501), the National Natural Science Foundation of China (40701032, 40801030), the National Essential Scientific Program of the Ministry of Science and Technology of China (2006FY110200), the Innovation Project of the Chinese Academy of Sciences (KZCX2-YWGJ04), a Grant-in-Aid for Science Research (19253001) from the Ministry of Education, Culture, Sports, Science and Technology of Japan, and Grants of Young Scientists A (20686033, 21686045) of the Japan Ministry of Education, Science and Culture and Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST). We also thank two anonymous reviewers and T. Scambos, Scientific Editor, for helpful comments and suggestions.
