Journal of Hydrology 444–445 (2012) 146–160

Contents lists available at SciVerse ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Catchment-scale reconstruction of glacier mass balance using observations
and global climate data: Case study of the Hailuogou catchment, south-eastern
Tibetan Plateau
Yong Zhang a,b,⇑, Yukiko Hirabayashi a, Shiyin Liu bab

Institute of Engineering Innovation, The University of Tokyo, 2-11-16 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan
State Key Laboratory of Cryospheric Sciences, Cold & Arid Regions Environmental & Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, Chinaa r t i c l e

i n f oArticle history:
Received 7 September 2011
Received in revised form 2 April 2012
Accepted 8 April 2012
Available online 17 April 2012
This manuscript was handled byKonstantine P. Georgakakos, Editor-in-Chief
Keywords:
Glacier
Mass balance
Debris cover
Catchment scale
Tibetan Plateaus u m m a r y
Debris-covered glaciers are common in the Tibetan Plateau, where ablation zones are mantled in a supraglacial debris cover that inﬂuences glacier mass balance, runoff, and response to climate change by affecting the melt rate of the underlying ice. The impact of debris cover has not yet been taken into account inregional- or global-scale assessments of glacier mass balances and freshwater resources by using physically based numerical models. Here, a surface energy–mass balance model that accounts for the signiﬁcance of debris cover and its effect on the ice melt rate is applied to reconstruct the glacier mass balanceof Hailuogou catchment, which is located in the south-eastern Tibetan Plateau and contains three debriscovered and four debris-free monsoonal maritime glaciers. According to our calculations, the glaciers in
Hailuogou catchment show a mean annual balance of 0.42 m water equivalent (w.e.) per year, for a total
mass loss of 24.3 m w.e. over the period 1952–2009. A comparison of summer temperature- and precipitation-mass balance/equilibrium line altitude (ELA) relations indicates that the glaciers in the catchment
are much more sensitive to temperature change than to precipitation change. In the last 20 yrs, increasing summer temperature is the main cause of rapid wasting of the glacier mass in the catchment.
Meanwhile, the presence of supraglacial debris markedly accelerates glacier mass loss, resulting in theunstable termini of debris-covered glaciers in Hailuogou catchment. This highlights the importance ofdebris cover for understanding glacier mass balance and hydrology in the Tibetan Plateau.
Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction
Mountain glaciers are a vital source of water in many of theworld’s river basins (Lemke et al., 2007) as well as sensitive indicators of climate change (Oerlemans and Fortuin, 1992). Glacier massbalance (mass gain or loss at the surface over a certain time period)
forms a vital link between climate change and glacier dynamics andhydrology, and its monitoring is the best way to infer climatechange from glaciers. The Tibetan Plateau and surrounding mountain ranges contain massive glaciers with a total area of
116,180 km2, whereas mass balance measurements have been
made on only 40 glaciers at one time or another since 1946 (Dyurgerov and Meier, 2005). Most regional mass balance estimates for the
Tibetan Plateau glaciers are based on extrapolation of only about 40
mostly short duration (<5 yrs in most cases) and unevenly distributed surface mass balance records (Dyurgerov and Meier, 2005;
Ohmura, 2006; Lemke et al., 2007) covering less than 1% of the total⇑ Corresponding author at: Institute of Engineering Innovation, The University ofTokyo, 2-11-16 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan. Tel.: +81 3 5841 0717.
E-mail address: zhyong@sogo.t.u-tokyo.ac.jp (Y. Zhang).
0022-1694/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.jhydrol.2012.04.014

glacier area. Furthermore, the Tibetan Plateau contains a large number of debris-covered glaciers. In these cases, ablation zones aremantled in a supraglacial debris cover that inﬂuences glacier melt
rates (Nakawo and Young, 1982; Mattson et al., 1993; Mihalcea
et al., 2006; Zhang et al., 2011) and terminus dynamics (Scherler
et al., 2011), which can modify a glacier’s response to climatechange. Although several numerical models have been proposedfor calculating the insulating effect of debris cover on different glaciers (e.g. Nakawo and Young, 1982; Fukushima et al., 1987; Braunet al., 1993; Nicholson and Benn, 2006; Reid and Brock, 2010; Zhang
et al., 2011; Immerzeel et al., 2012), so far no physically based
numerical mass-balance models have been developed on a large
spatial unit (e.g. on a catchment scale) that can consider the significance of debris cover and its inﬂuence on ice melt rates that havedetermined glacier mass balances. We argue that the estimates reported in previous studies (e.g. Dyurgerov and Meier, 2005; Ohmura, 2006; Lemke et al., 2007) lead to substantial uncertainty inthe mass loss from the Tibetan Plateau, especially in regions wheredebris-covered glaciers exist. Hence, there is an urgent need for better methods, which can consider the signiﬁcance of debris cover andits impact on the ice melt rate, to quantify mass changes from the
Tibetan Plateau glaciers.

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

In recent years, physically based mass balance models havebeen applied to individual glaciers driven using in situ meteorological observations (e.g. Arnold et al., 1996; Fujita and Ageta, 2000;
Klok and Oerlemans, 2002; Hock and Holmgren, 2005; Sakai
et al., 2009). However, the application of these models to large spatial units (e.g. on a catchment scale) has been limited because of
the inadequate spatial and temporal coverage of local meteorological records and mass balance observations (Barry, 2006). To overcome this problem, climate model simulations have been used. A
few studies have investigated the ability of climate model simulations to drive energy balance models (e.g. Reichert et al., 2001;
Bougamont et al., 2005; Hock et al., 2007; Fujita et al., 2011; Fujita
and Nuimura, 2011), and degree-day models (e.g. Hock et al., 2007;
Hirabayashi et al., 2010) for modeling mass balance. Their results
indicate that climate model simulations provide an opportunity
to be used for reconstructing mass balance over a long period in
the regions that are not routinely monitored.
In this study, we implement a surface energy–mass balance modeling approach that takes into account the signiﬁcance of debris cover and its impact on the ice melt rate, to reconstruct the glacier massbalance of Hailuogou catchment in the south-eastern Tibetan Plateau, forcing the model by observations and global climate data.
Hailuogou catchment is chosen as the best case because it offers an
opportunity to study a monsoonal maritime glacier system withdebris-covered and debris-free glaciers in the south-eastern Tibetan
Plateau, where speciﬁc, though incomplete, information is available
about both glaciology and meteorology. The main aim of this study isto present a simple, physically-based approach considering the signiﬁcance of debris cover and its inﬂuence to estimate as accurately as
possible the long-term variation and trend in surface mass balanceon a catchment scale. By this approach, the effect of supraglacial debris cover on glacier mass balance can be assessed. This, in turn, allows
us to document ongoing glacier and climate changes in the southeastern Tibetan Plateau, one of the largest and most sensitive
remaining areas of temperate ice (Li and Su, 1996), covering approximately the same area as the Southern Patagonia Iceﬁeld (Rignot
et al., 2003). Speciﬁcally, we validate the global climate data in the
Hailuogou catchment and explore the potential to use global climate
products for mass balance modeling of glaciers on the Tibetan Plateau, where the harsh climatic conditions and remoteness hamper
ground-based monitoring and result in very poor data coverage.

2. Study area
Hailuogou catchment (80.5 km2) is located on the eastern slopeof Mount Gongga, south-eastern Tibetan Plateau (29.6°N, 101.9°E)(Fig. 1). The catchment contains seven glaciers covering an area of
36.44 km2 (Table 1) (Pu, 1994). These seven glaciers are named from
Glacier No. 1 to Glacier No. 7 in Fig. 1. Hereafter, the largest glacier(No. 3) is called Hailuogou Glacier. Among the seven glaciers, theablation zones of Glacier No. 1, No. 2, and Hailuogou Glacier are covered by supraglacial debris cover supplied from mixed ice/snow/
rock avalanches occurring on the icefall and the surrounding walls
through frost weathering processes and structural rockfalls (Li andSu, 1996). In situ debris thickness surveys on Hailuogou Glacierindicated that the debris thickness increases from just several
millimetres of patchy cover below 3650 m above sea level (a.s.l.)to >1 m at the glacier terminus (Zhang et al., 2011).
The catchment climate is characterized by the southeast monsoon in summer and westerly circulation in winter (Li and Su,
1996). The mean annual air temperature (1988–2004) observedat the Gongga Alpine Ecosystem Observation and Research Station(hereafter, GAEORS; 3000 m a.s.l.; Fig. 1) is 4.1 °C, and mean annual
precipitation (1988–2004) is 1906 mm w.e. Summer precipitation(May–October) accounts for 80% of the annual total.

147

Glaciers in the catchment have retreated considerably since theearly 20th century (Table 1) (Su et al., 1992; Liu et al., 2010; Zhang
et al., 2010). The total glacier area of the catchment decreased by
2.9 km2 over the period 1966–2007, representing 7.8% of the total
glacier area in 1966 (Liu et al., 2010). Zhang et al. (2010) noted thatthe ablation area of Hailuogou Glacier thinned by 1.1 ± 0.4 m yr1during the period 1966–2009, and the thinning has accelerated signiﬁcantly since 1989. The termini of Glacier No. 2 and Hailuogou
Glacier retreated by 1100 m during the period 1930–1990 (Suet al., 1992) and 1878 m during the period 1930–2007 (Suet al., 1992; Liu et al., 2010), respectively.
3. Data
We use various datasets, including meteorological observationsat GAEORS, daily global gridded data of temperature and precipitation from the closest grid of the catchment, runoff records fromGAEORS, ablation data from Hailuogou Glacier, glacier area data,
debris thickness data, and a digital elevation model (DEM) with aresolution of 30 m. These datasets are brieﬂy described below.
3.1. Meteorological data
GAEORS (Fig. 1) has been operated near the terminus of
Hailuogou Glacier since 1988. Daily temperature, wind speed, relative humidity, and precipitation data are available at GAEORS forthe period 1988–2009. An automatic weather station (AWS; Fig. 1)
was installed in the ablation area of Hailuogou Glacier at 3550 ma.s.l. in April 2008. Air temperature, relative humidity, wind speed,
and the incoming and reﬂected solar radiation are recorded at 10min intervals. Precipitation has been measured by a precipitation
gauge installed at the margin of Hailuogou Glacier (3500 m a.s.l.;Fig. 1) since April 2008, and its 30-min average was stored. All variables were compiled to daily values.
To reconstruct the long-term glacier mass balance variation in
Hailuogou catchment, a global 0.5° gridded dataset of near-surface
temperature and daily precipitation for 57 yrs from 1951 to 2007 isused. The temperature dataset was generated by Hirabayashi et al.
(2008) based on monthly temperature values of the Climate
Research Unit version TS 2.1 data set for the period 1951–2002(CRU; Mitchell and Jones, 2005), and the dataset of Fan and van
den Dool (2008) for the period 2003–2007 (hereafter, H08). Full details of the H08 dataset and generation methods have been reported
by Hirabayashi et al. (2008). The daily precipitation dataset for
1951–2007 was created by collecting rain gauge observation data
across Asia (Yatagai et al., 2009). This is the only long-term continental-scale daily product that contains a dense network of dailyrain gauge data for Asia including the Himalayas and mountainous
areas in the Middle East (Yatagai et al., 2009) (hereafter, Y09). Fulldetails of the Y09 dataset can be found in the work of Yatagai
et al. (2009) or at http://www.chikyu.ac.jp/precip/index.html. The
altitude information of each grid cell is derived from an SRTM digital
elevation model (February 2000; Jarvis et al., 2008).
3.2. Runoff and area data
Runoff derived from the catchment has been measured at theriver outlet of the catchment located at the terminus of Hailuogou
Glacier since 1994 (Fig. 1). Records are available for the period
1994–2002. In this study, monthly runoff observations are usedas simulation targets to calibrate and verify the energy–mass balance model.
The area of each glacier in the catchment derived from the Chinese Glacier Inventory (Pu, 1994) is shown in Table 1, representing
the area information before 1966. To estimate the glacier area in148

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Fig. 1. Map of Hailuogou catchment on the eastern side of Mount Gongga in the south-eastern Tibetan Plateau, spatial distribution of thermal resistance of the debris layer
and area-altitude distributions of the catchment, glacier, and debris-covered surface. GAEORS denotes the Gongga Alpine Ecosystem Observation and Research Station(square). AWS denotes the automatic weather station (star). The circle marks the hydrological station for runoff observation. The cross is the precipitation gauge. Glacierboundaries in 1966 and 2007 are derived from Liu et al. (2010). Spatial distribution of thermal resistance of the debris layer is derived from Zhang et al. (2011). The histogramshows the area-altitude distributions of the debris-covered ice (gray), the glaciers (red), and the catchment (blue). (For interpretation of the references to color in this ﬁgurelegend, the reader is referred to the web version of this article.).

Table 1
Characteristics of seven glaciers in Hailuogou catchment. Glacier area, length and elevation are derived from the Chinese Glacier Inventory (Pu, 1994). Annual variation rate (AVR)of glacier area is derived from Liu et al. (2010).
Glacier

Elevation (m a.s.l.)Length (km)Area (km2)

AVR (% yr1)

Debris-covered

No. 1
No. 2
No. 3
No. 4
No. 5
No. 6
No. 7

3740–5400
3600–6368
2990–7556
5120–5840
5080–6000
4720–6136
5280–5800

4.3
5.6
13.1
0.9
1.1
2.0
0.8

1.44
7.67
25.71
0.33
0.37
0.66
0.26

0.38
0.21
0.09
1.80
1.68
1.45
1.52

Yes
Yes
Yes
No
No
No
No

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

different periods, glacier outlines in 1966, 1975, 1994, and 2007
are used. These were created by Liu et al. (2010) based on topographic maps and remote sensing data. Full details of these data
can be found in the work of Liu et al. (2010).
3.3. Mass balance data
The ﬁrst ablation measurement in the catchment was performed on the ablation area of Hailuogou Glacier in 1982 with23 stakes (Li and Su, 1996). Since then, ablation measurements
have been conducted discontinuously at a network of ablation
stakes distributed over the ablation area of Hailuogou Glacier(Liu et al., 1994; Su et al., 1996; Zhang et al., 2011). Ablation data
in 1982/83 (Li and Su, 1996), 1990/91–1993/94 (Liu et al., 1994; Su
et al., 1996), and 2008 (Zhang et al., 2011) are used to verify theperformance of the energy–mass balance model.
Snow depth was measured at 150 points across an altitudinalrange of 4900–5500 m a.s.l. in October 1990 (Aizen and Aizen,
1994; Liu et al., 1994). A mean snow density of 415 kg m3 was obtained from ﬁve snow pits and a crevasse with distinct annual layer
in the accumulation area of Hailuogou Glacier (Aizen and Aizen,
1994; Liu et al., 1994). This snow density is used to convert allsnow-related mass balance measurements into water equivalent.
3.4. Debris thickness data
The ﬁrst debris thickness surveys in the catchment were carriedout on Hailuogou Glacier in 1982 at 23 points across an altitudinalrange of 2960–3600 m a.s.l. (Li and Su, 1996). In 2009, extensive
in situ surveys of debris thickness were conducted at 300 points
on Hailuogou Glacier, and these points were surveyed across theentire ablation area in the altitudinal range of 2900–3600 m a.s.l.
(Zhang et al., 2011).
Data of thermal resistances of debris layers are used to map the
spatial distribution of debris cover in the catchment (Fig. 1) andcalculate the ice melt rate beneath the debris. The thermal resistance of the debris layer is deﬁned as debris thickness divided by
the thermal conductivity of the debris layer (Nakawo and Young,
1981, 1982). Field determinations of the thickness and thermal
conductivity of the debris layer are especially time-consuming
and not practical on glaciers. Zhang et al. (2011) calculated thermal
resistance of debris layers from the surface temperature derived
from the thermal infrared bands of ASTER and the net radiation
estimated from the visible/near-infrared bands of ASTER and
NCEP/NCAR reanalysis data that corresponded to the nearest timeand location of ASTER acquisition. An analysis of the transverse andlongitudinal proﬁles of the ground-surveyed debris thicknesses
and ASTER-derived thermal resistances of debris layers indicatedthat the ground-surveyed debris thicknesses and ASTER-derived
thermal resistances of debris layers correlate strongly over the entire ablation area, and across- and along-glacier patterns of ASTERderived thermal resistance correspond well with spatial patterns ofdebris thickness, which can reﬂect large-scale variations in the extent and thickness of the debris cover (Zhang et al., 2011). Full details of the thermal resistance of the debris cover and its
calculation method can be found in the work of Zhang et al. (2011).
3.5. Digital elevation model
A DEM of the catchment was produced from aerial photographs
acquired in 1989 (Zhang et al., 2010). The accuracy of the DEM wasvalidated by the global positioning system (GPS) data on glacierfree terrain, yielding a root mean square error (RMSE) of 11.2 m(Zhang et al., 2010). We use the DEM to spatially distribute meteorological data and compute the area of each elevation band in thecatchment.

149

4. Methods
4.1. Catchment discretization
Hailuogou catchment is divided, at intervals of 50 m, into a setof elevation bands based on the DEM. Each elevation band is assumed to have a homogenous hydrological behavior. The areas ofeach elevation band in 1966, 1975, 1994, and 2007, which are used
to calculate the area-averaged mass balance, are obtained from theDEM based on glacier outlines for 1966, 1975, 1994, and 2007. The
area-altitude distributions of the glacier and the catchment are
shown in Fig. 1. Although shrinkage of the glaciers was reported
in the catchment (e.g. Su et al., 1992; Liu et al., 2010; Zhang
et al., 2010), glacier altitude-area distributions in 1966, 1975,
1994, and 2007 are correspondingly supposed to be constant forthe periods 1952–1966, 1967–1975, 1976–1994 and 1995–2009.
Ablation zones of Glacier No. 1, Glacier No. 2, and Hailuogou
Glacier are subdivided into debris-covered and debris-free surfaces. This separation is completely based on the spatial distribution of the thermal resistance of the debris layer (Fig. 1). The
thermal resistance of the debris layer can be used to be as a proxy
for the debris thickness over large areas because the thermal resistance of the debris layer correlates strongly with the ground-surveyed debris thickness, the spatial pattern of which reﬂects
large-scale variations in the extent and thickness of the debris cover (Zhang et al., 2011). The area-altitude distribution of debris-covered surfaces is shown in Fig. 1.
4.2. Meteorological data pre-processing
For each elevation band, the temperature and precipitation timeseries are interpolated according to the mean elevation. The airtemperature is assumed to decrease with increasing altitude witha constant lapse rate. Cao and Cheng (1994) estimated the temperature lapse rate at Mount Gongga as 4.3 °C km1 on the eastern
slope and 5.4 °C km1 on the western slope. Here, the lapse rate
is ﬁxed to 4.3 °C km1. This rate is conﬁrmed by recent observations at GAEORS and the AWS in the ablation area of Hailuogou
Glacier.
Precipitation increases linearly with altitude in this region(Aizen and Aizen, 1994; Cao and Cheng, 1994; Cheng, 1996; Liuet al., 2010). We apply a precipitation gradient dprec (% of precipitation increase per meter of elevation increase) from the snout to the
top of the glacier. Thus, the precipitation, P, at each elevation bandis calculated as

P ¼ PGAEORS ½1 þ dprec ðh  hGAEORS Þ;ð1Þ

where PGAEORS is precipitation at GAEORS, and h and hGAEORS are thealtitudes of the elevation band and GAEORS, respectively. Based on
results of previous studies (Aizen and Aizen, 1994; Cao and Cheng,
1994; Liu et al., 1994, 2010; Cheng, 1996), we estimate the daily
orographic precipitation gradients of this region varying between
12% and 60% per 100 m. However, precipitation has spatial variability in the glacierized regions of high mountainous catchments.
Accurate spatial interpolation of precipitation data series is therefore difﬁcult to obtain and the prevailing extreme conditions imply
signiﬁcant measurement uncertainty. In particular, the value of the
precipitation gradient is difﬁcult to estimate in the accumulation
zones of the catchment, where precipitation data are scarce. This
difﬁculty in the spatial interpolation of precipitation is a signiﬁcantsource of modeling uncertainty for glacier mass balance and runoff
simulations. Considering this, here we treat the precipitation gradient as a tunable model parameter.
In addition, because wind speed only weakly correlates with themelt rate (Ohmura, 2001), and has a little inﬂuence on the mass

150

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

balance results (Fujita and Ageta, 2000; Sakai et al., 2009; Fujita
et al., 2011), wind speed observed from GAEORS is not interpolated
and assumed representative for the entire area calculated.
4.3. Global climate data validation and bias correction
Air temperature and precipitation observed from GAEORS areused to validate the global climate data (H08 and Y09). A linear
regression analysis is used to determine the correlation between
the global climate data and observed variables on daily, monthly,
seasonal and annual timescales. The H08 temperature is validated
for the period 1988–2007 using data from GAEORS, and the Y09
precipitation is validated for the period 1988–2004. After evaluating the association between the two datasets, we correct the global
climate data for model bias on a daily timescale by a linear regression method because the global climate data are frequently biased
relative to in situ observations due to their parameterization
schemes and coarse grid resolution. Linear regression equations
are established between in situ meteorological observations andglobal climate data over the period 1988–1997. These are expressed as follows:

T GAEORS ¼ 0:995T H08  2:559 ðr ¼ 0:88; p < 0:001;
RMSE ¼ 2:7  CÞ;

ð2Þ

PGAEORS ¼ 1:250PY09 þ 2:199 ðr ¼ 0:71; p < 0:001;
RMSE ¼ 0:003 mÞ;ð3Þwhere TGAEORS and PGAEORS are daily temperature and precipitation
at GAEORS, respectively, and TH08 and PY09 are H08 temperature
and Y09 precipitation, respectively.
4.4. Surface energy–mass balance model
The surface energy–mass balance model consists of two coupled
components. The ﬁrst is to compute the energy available for melting
from energy exchange between the debris-covered/free surface andthe atmosphere, and accumulation. The second is to treat processes
occurring in the subsurface after meltwater percolates in the underlying layers (Fig. 2). The model calculates the surface mass balance
for each elevation band at a time step of one day. The annual massbalance for each elevation band (i = 1. . .n), bi, is calculated as

bi ¼

t1
X
½CaðtÞ þ MðtÞ þ RF ðtÞ;

ð4Þ

t¼t 0

where M represents ablation (mass loss is deﬁned negatively), Ca(positive sign) is accumulation (snow), and RF (positive sign) is
refreezing. All terms are in units of m w.e. t0 is the start date of
the mass balance year (here, 1 October), and t1 is the end date of
the mass balance year (30 September in the following year). The annual mass balance of the entire catchment, B, is estimated as thearea-weighted sum of the mass balance in all elevation bands,
which is calculated as

Pn
bi Si
B ¼ Pi¼1;n
i¼1 Sið5Þwhere Si is the discrete value of the surface area for each elevation
band (i = 1. . .n) (m2).
The glacier melt, M, in the ith elevation band is computed using
a surface energy-balance model. For the debris-free ice, the energy
available for melting, QM, is calculated fromQ M ¼ ð1  aÞRS þ RLd þ RLu þ Q S þ Q L þ Q R þ Q G ;ð6Þ

where RS and RLd are the downward short- and long-wave radiation
ﬂuxes, respectively, RLu is the upward long-wave radiation ﬂux, QS,

QL, QR and QG are the net sensible heat ﬂux, net latent heat ﬂux, heatﬂux by rain on the surface and conductive heat ﬂux into the glaciersurface, respectively, and a is the surface albedo. All terms are taken
to be positive toward the surface, in units of W m2. Zhang et al.
(2011) analyzed the effect of QR on ice melt rates on Hailuogou Glacier, and found that it was not signiﬁcant. Hence, the heat ﬂux by rain
is not taken into account. Downward long-wave radiation is obtained
from the air temperature, relative humidity and ratio of downward
short-wave radiation to that at the top of the atmosphere using anempirical scheme (Fujita and Ageta, 2000). The upward long-wave
radiation is calculated from the Stefan–Boltzmann constant andthe surface temperature assuming a black body for the snow/ice surface, and turbulent sensible- and latent-heat ﬂuxes were calculatedusing the bulk method. They are expressed as follows

RLu ¼ erðT S þ 273:2Þ4
Q S ¼ ca qa CUðT a  T S Þ

ð7Þ

Q L ¼ le qa CU½rhqðT a Þ  qðT S Þ;where r is the Stefan–Boltzmann constant (5.67  108 W m2 K4),
e is the emissivity of the glacier surface (taken to be 1), TS is the surface temperature (°C), Ta is the air temperature (°C), ca is the speciﬁcheat capacity of air (1006 J kg1 K1), qa is the density of air (kg m3),
C is the bulk coefﬁcient for sensible and latent heat (0.002), U is thewind speed (m s1), le is the latent heat of the evaporation of water(2.5  106 J kg1), q is the saturated speciﬁc humidity (kg kg1) calculated as a function of air temperature using the empirical equations presented by Kondo (1994), and rh is the relative humidity.
Before 1988, relative humidity data were not available in the catchment. Fig. 3a shows the relation between the monthly means of precipitation and relative humidity established based on observations at
the AWS during the period October 2008–July 2009. The relation
indicates relative humidity increasing with precipitation and is used
to estimate the daily relative humidity from the daily precipitation.
As shown in Fig. 3b, the estimated relative humidity is in agreementwith the observed data for the period May–September 2008.
The absorbed short-wave radiation is calculated from thesurface albedo and downward short-wave radiation. The surface
albedo calculation in the model follows the albedo scheme of Fujita(2007). This method was applied to different regions and validated
against observations (e.g. Sakai et al., 2009; Fujita et al., 2011;
Fujita and Nuimura, 2011; Yasunari et al., 2011). In their studies,
the albedo is calculated by a simpliﬁed concept of multiple scattering in the ice plate whose thickness is related to the surface snow
density. The albedo, a, in the Fujita (2007) scheme is given by

a ¼ rI þ

ð1  r I Þ2 s;
1  rI s

ð8Þ

where

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ
22

I Þ RI
s ¼ ð1T I Þ Rð1TI
2

I Þ expðkI lI Þ
T I ¼ ð1r
1r 2 expð2k l Þ
I

I I

2
I Þ r I expð2kI lI Þ
RI ¼ r I þ ð1r
1r 2I expð2kI lI Þ

In these equations rI is the reﬂectivity of ice (0.018), RI is thereﬂectance of a single ice plate, TI is the transparency of a single
ice plate, lI is the effective ice thickness, and kI is the absorption
coefﬁcient of ice (assumed to be 10 m1). The effective ice thickness lI can be determined from

lI ¼

2

qi S

;

ð9Þ

where qi is the density of ice (assumed to be 900 kg m3), and S⁄(m2 kg1) represents speciﬁc surface area (the area of the surfaceof the ice particles in unit mass of snow), computed with

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

151

Fig. 2. Schematic structure of the surface energy–mass balance model. ST denotes the surface temperature. DSEB and DFSEB denote the surface energy balance of debriscovered and debris-free surfaces, respectively.

log10 S ¼ 15:32  109 q3 þ 16:65  106 q2  7:30
3

 10 q þ 3:23;

ð10Þ

In Eq. (10), a meltwater effect is invoked when the snow temperature rises to 0 °C, and S⁄ will be decreased to 60% of its originalvalue. The surface snow density q changes with snow densiﬁcation, which is estimated using a model of snow densiﬁcation established by Motoyama (1990). He suggested that the snow layer isdensiﬁed due to viscous compression. The temporal variation insnow density is calculated from a compactive viscosity factor estimated by an empirical formula with respect to the density of asnow layer at a given depth and the overburden load exerted onthat layer. The inﬂuence of water, which enhances snow compaction, is taken into account by changing the compactive viscosity
factor for wet snow. The albedo of debris-free ice was set to 0.3(Zhang et al., 2011), whereas that of the new snow was assumed
to be 0.8 (Li and Su, 1996).
The downward short-wave radiation is calculated by using amethod developed by Matsuda et al. (2006). They found a favorable correlation between precipitation and atmospheric transmissivity of solar radiation in terms of the monthly values on
different glaciers of the Tibetan Plateau. This method was applied
to the July 1st Glacier on the Tibetan Plateau where it performed
well (Sakai et al., 2009). The relation between the monthly meantransmissivity and the monthly mean precipitation observed at

the AWS during the period October 2008–July 2009 is shown in
Fig. 3c. This is used to estimate the daily mean transmissivity, s,from daily precipitation data. The downward short-wave radiationis then obtained from the daily mean transmissivity and theincoming short-wave radiation at the top of the atmosphere. As
shown in Fig. 3d, short-wave radiation ﬂuctuations are well simulated for the period May–September 2008, although the estimatedshort-wave radiation is greater than the observations that are<110 W m2. This is because the relation is established based onmonthly data, for which low transmissivity data with high precipitation are not taken into account. In addition, glacier melt, whichmainly determines the glacier mass balance, occurs when the
short-wave radiation is large. Thus, overestimation of low
short-wave radiation has little effect on the glacier mass balance
estimate.
The surface temperature, TS, is obtained to satisfy all heat balance equations by iterative calculation of the conductive heat ﬂux,
which is obtained by calculating the temperature proﬁle of thesnow layer and/or glacier ice. The melt rate of the debris-free ice,
M, is calculated as

M¼

QM

qi Lf

;

where Lf is the latent heat of fusion (3.34  105 J kg1).

ð11Þ

152

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Fig. 3. (a) Relation between monthly precipitation and relative humidity observed from October 2008 to July 2009, (b) relative humidity observed from May to September
2008 and estimated using the relation shown in (a), (c) relation between monthly precipitation and monthly average transmissivity of solar radiation (downward short-wave
radiation/short-wave radiation at top of the atmosphere), s, observed from October 2008 to July 2009, and (d) daily solar radiation observed from May to September 2008 andits estimation using the relation shown in (c).

The catchment contains three debris-covered glaciers. The total
debris-covered area accounts for 39% of the total ablation zone
area of three glaciers. Fig. 1 shows the altitude-area distributionof debris-covered surfaces. At the debris-covered surface, the energy balance is expressed as

Q 0G ¼ ð1  a0 ÞRS þ RLd þ RLu þ Q S þ Q L þ Q R ;ð12Þ

where Q 0G is the conductive heat ﬂux into the debris layer, and theother energy components and calculations are the same as thosementioned above. All terms are taken to be positive toward the debris surface, in units of W m2. The albedo of the debris-covered surface, a0 , was set to 0.2 (Zhang et al., 2011). The only heat ﬂux
considered to reach the glacier ice is the conductive heat ﬂux withsimplifying assumption of a linear temperature proﬁle within thedebris layer and the constant heat ﬂux stored in the debris layerfrom day to day. Therefore, the heat ﬂux used for the ice melt,
QM, can be calculated as

Q M ¼ Q 0G ¼

ðT 0S  T I Þ
;
R

ð13Þ

where R is the thermal resistance of the debris layer (m2 K W1), T 0S
is the debris surface temperature (°C), and TI is the ice temperature
at the ice-debris surface, which was set to 0 °C (Zhang et al., 2011).
All energy-balance components at the debris-covered surface,
except the short-wave radiation term, are explicitly determinedfrom the debris surface temperature. The debris surface temperature, T 0S , can be calculated numerically by iteration. The resulting
daily mean surface temperature is then used to calculate the rateof ice melt from Eqs. (11) and (13) for the prescribed thermal resistance of the debris layer. Zhang et al. (2011) found that the thermal
resistance of the debris layer at a speciﬁc site could be regarded as
constant usually in a melting season on Hailuogou Glacier.

Therefore, the thermal resistance of the debris layer was assumed
to be constant over the study period. Full details of the entire model have been reported by Zhang et al. (2011).
Some of the meltwater can be ﬁxed to the glacier by refreezing if
the glacier ice is cold enough (Fujita et al., 1996). Our calculations of
the refreezing amount follow Fujita and Ageta (2000), which calculate the refreezing amount from considering the conduction of heat
into the glacier ice and the snow layer and the presence of water at
the interface between the snow layer and glacier ice. The refreezing
amount RF in the Fujita and Ageta (2000) scheme is given by

RF ¼

qi ci
Lf

Z zc
interface

DT Z dz þ

qs c iLf

Z interface

DT Z dz;

ð14Þ

surfacewhere qs is the average density of snow (415 kg m3) (Aizen and Aizen, 1994; Liu et al., 1994), ci = 2100 J kg1 K1 is the speciﬁc heat
capacity of ice, z is the depth from the surface or snow-ice interface(m), zc is the depth where the annual amplitude of ice temperatureis <0.1 °C (assumed to be 20 m from the snow-ice interface), and
DTZ is the temperature difference of the snow/ice at a given depth
during a given period (°C). In addition, refreezing during winterand shorter cooling events is considered by assuming the maximum
water content as 5% in volume.
For each elevation band, snow accumulation is modeled fromthe precipitation value by using a simple temperature thresholdto discriminate rain from snow. The threshold temperature was
set to 2.0 °C (Liu et al., 2009). A mixture of snow and rain is assumed for a transition zone ranging from 1 K above and 1 K below
the threshold temperature. Within this temperature range, thesnow and rain percentages of total precipitation are obtained fromlinear interpolation. Redistribution of the original snowfall by windtransport or avalanches is not taken into account.

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

4.5. Model calibration

5. Results

Model calibration involves constraining model parameters to get
the best ﬁt between the model simulations and the observed data.
The parameters are adjusted within the accepted ranges according
to one or more criteria that assess the goodness-of-ﬁt between thesimulated results and the observed data. In this study we use the observed runoff to optimize the precipitation gradient within therange from 12% to 60%. The calibration is performed to minimize
the difference between observed and modeled runoff in different
years. The precipitation gradient associated with the best agreement between observed and modeled runoff in different years is
chosen. The presented model enables the estimation of the runoff
from the catchment. The runoff contributions of all elevation bands
are added to the total discharge at the outlet of the entire catchment. No routing between the spatial units and the river outlet is
performed because the catchment has rather steep slopes (Fig. 1).
This means that the runoff delay due to routing in the river networkis thus much smaller than the given time step of one day. Liu et al.
(2010) analyzed the water balance of the catchment using a hydrological method, and their results indicated that their method for calculating water balance of the catchment performed reasonably wellwithout considering groundwater system and changes in evaporation and/or sublimation over the glacier surface, and the main contributions to the runoff were meltwater and rainfall. Furthermore,
the runoff in the catchment shows a continuously increasing trend
for the period 1994–2002, especially after 1997. The annual runoff
reaches the maximum in 2002, which increases by 42% comparedto the mean annual runoff over the period 1994–2002. Therefore,
our model is calibrated against observed runoff in 1997 and 2002.
The runoff, Q (m3 s1), can be calculated as

5.1. Validation of global climate data

Q ¼ Q1 þ Q2
Q 1 ¼ M þ Pliq  RF ;

ð15Þ

Q 2 ¼ P tot  E
where Q1 and Q2 are the runoff from the glacier and glacier-freezone, respectively, Pliq is rain on the glacier, Ptot is the total precipitation, and E is the evaporation in the glacier-free zone calculated
from the temperature and precipitation using the empirical equation of Kang et al. (1999). In addition, we assume there is no
remaining snow at the end of the melt season, since our aim is tosimulate the monthly runoff. Therefore, precipitation is not separated into rain and snow in the glacier-free zone.
Most mass balance models presented in previous studies arecalibrated to get a good ﬁt between observed and simulated variables by attempting to minimize a single measure of model performance (e.g. Oerlemans and Fortuin, 1992; Braithwaite and Zhang,1999; Hock et al., 2007). However, there is no best criterion for statistically evaluating the quality of models because all models are
simpliﬁcations of reality. Hence, calibration using a single criterion
often constrains the model parameters to ﬁt certain characteristicsof the observed data, but neglects the remaining features. Here, we
adopt multiple assessment criteria to assess the model performance. With the exception of the correlation coefﬁcient r, threeassessment criteria are used: the efﬁciency criterion R2 (Nashand Sutcliffe, 1970), a modiﬁed version R2ln (Hock, 1999), and theroot mean square error (RMSE). The criteria R2 and RMSE havebeen widely used to measure the quality of ﬁt between simulations
and observations. The criterion R2ln is a measure of the agreement
between simulations and observations emphasizing low-ﬂow conditions. When the logarithms are found, discrepancies between
modeled and observed runoff during peak ﬂows exert less impact
on the numerical value. The criteria, R2 and R2ln , assume a valueof 1 in case of perfect agreement.

153

Unless otherwise stated, all signiﬁcance levels of correlation
coefﬁcients (r) presented herein are p < 0.001. As shown in Table
2, daily, monthly, and annual H08 temperature and Y09 precipitation means correlate well with the corresponding observations
from GAEORS, although correlation coefﬁcients (r) for precipitation
are lower than those for temperature. Daily, monthly, and annual
H08 temperature means correlate well with the observations fromGAEORS for the period 1988–2007, yielding coefﬁcients of 0.90,
0.98, and 0.71, respectively, while analysis of the daily, monthly,
and annual Y09 precipitation data yields correlation coefﬁcients(r) of 0.70, 0.96, and 0.53, respectively. Seasonally, the correlations
between the H08 and observed temperature are higher during thespring (March, April, May (MAM)), summer (June, July, August(JJA)), and autumn (September, October, November (SON)), compared to the winter (December, January, February (DJF)) (Table2). Further analysis of the results reveals that the winter correlation is unduly inﬂuenced by the Decembers of 5 yrs (1998–2002)
for which daily H08 temperatures are signiﬁcantly larger thanthe observations, while the H08 and observed temperatures in
December of other years show similar ﬂuctuations, although there
is a 433 m difference in elevation between the grid cell and
GAEORS. If we remove the inﬂuence of the ﬁve exceptional Decembers, the correlation coefﬁcient r during the winter rises to 0.85.
This indicates that the H08 simulates winter temperatures reasonably well for the majority of years, and is able to capture theobserved temporal variability well. Meanwhile, the Y09 precipitation data better capture the observed variability in spring (0.90),
summer (0.70) and autumn (0.95) than in winter (0.53) (Table 2).
Precipitation from spring to autumn accounts for 96% of the annual
total in Hailuogou catchment, based on the observation during
1988–2004. The glaciers in the catchment are the summer-accumulation type glaciers (Shi and Liu, 2000). Thus, despite the lower
correlation coefﬁcient (r) during the winter, precipitation in wintermay have little inﬂuence on the model performance.
Global climate data have to be corrected before use at the subgrid scale in the mass balance models. Hence the global climate
data series used in this study are bias-corrected by Eqs. (2) and(3) which are established between in situ meteorological observations observed data and the global climate data over the period
1988–1997. Fig. 4a and b shows the bias-corrected and observed
daily temperature averaged over the period 1998–2007 and precipitation averaged over the period 1998–2004. The correlation
coefﬁcient (r) between the bias-corrected and observed daily temperature is 0.96 with the RMSE value of 1.54 °C, and the correlation
coefﬁcient (r) between the bias-corrected and observed daily precipitation is 0.87 with the RMSE value of 0.001 m. Fig. 4c and dshow the variations in the bias-corrected and observed monthly
temperature over the period 1998–2007 and annual precipitation
over the period 1988–2004. As can be seen, bias-corrected global
climate data yields good agreement with the observations.
Overall, the above analysis suggests that the global gridded
temperature and precipitation data used in this study correspond
sufﬁciently well with the observations at GAEORS to be used as input for mass balance model in Hailuogou catchment.

5.2. Model calibration and veriﬁcation
Fig. 5 shows the best solution found for the runoff simulation indifferent years, and the ﬁt between observed and modeled runoff,
as measured by the efﬁciency criterion R2 , its modiﬁed version R2ln
and the RMSE. The results of model calibration show that there is

154

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Table 2
Correlation coefﬁcients (r) between the global climate data and observations fromGAEORS. All signiﬁcance levels of correlation coefﬁcients presented herein are
p < 0.001.

Daily
Monthly
AnnualSpring (May, April, May (MAM))
Summer (June, July, August (JJA))
Autumn (September, October, November(SON))
Winter (December, January, February (DJF))

Temperaturea

Precipitationb

0.90
0.98
0.71
0.98
0.91
0.99

0.70
0.96
0.53
0.90
0.70
0.95

0.64

0.53

a

Data were used for the period 1988–2007.
b
Data were used for the period 1988–2004.

close agreement between observed and modeled monthly runoff,yielding R2 > 0.9, R2ln > 0.7 for using observed meteorological data,
and R2 > 0.9, R2ln > 0.9 for using global climate data (Fig. 5). The runoff simulations in 1997 and 2002 by using observed meteorological
data and global climate data have small RMSE values compared tothe corresponding mean monthly discharge of 8.1 m3 s1 and13.7 m3 s1 (Fig. 5). This indicates that the precipitation gradientof 18% per 100 m can give a good agreement between observed

and modeled runoff, and is therefore selected as the best solution
for further analysis.
The calibrated model (Fig. 5) is veriﬁed against an independent
set of runoff data collected over the periods 1994–1996 and 1998–
2001. Scatter diagram of observed and modeled monthly runoff isshown in Fig. 6, along with a plot of the corresponding cumulative
variation. Overall, the calibrated model performs well in the validation period. The efﬁciency criterion R2 and its modiﬁed version
R2ln are respectively 0.84 and 0.80 for observed meteorological data,
and 0.84 and 0.86 for global climate data (Fig. 6a). The monthly
runoff simulations using observed meteorological data and global
climate data have small RMSE values compared to the mean
monthly discharge of 13.4 m3 s1, which equal to 2.68 and2.7 m3 s1, respectively. As shown in Fig. 6b, the model performs
well with respect to the cumulative runoff over the period 1994–
2002. The cumulative runoff simulated using observed meteorological data is underestimated by 10% over the period
1994–2002, whereas that simulated using global climate data isoverestimated by 4% at the same period. As a result, the modeled
cumulative runoff can capture the long-term trend.
The calibrated model allows for the effect of debris cover on glacier ablation, which is veriﬁed by comparing the model results toobservations measured at stakes on Hailuogou Glacier in different
periods. Scatter diagrams of modeled and measured ablation at

Fig. 4. (a) Observed daily air temperature and bias corrected H08 temperature averaged over the period 1998–2007, (b) observed daily precipitation and bias-corrected Y09
precipitation averaged over the period 1998–2004, (c) comparison of observed monthly temperature and bias-corrected H08 temperature over the period 1998–2007, and (d)comparison of observed and corrected annual precipitation over the period 1998–2004.

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

155

Fig. 5. Observed and modeled runoff of the catchment in 1997 (a) and 2002 (b) using the optimized precipitation gradient with respect to the assessment criteria (R2, R2ln and
RMSE (m3 s1)). Qmod, observed data and Qmod, global data denote the simulations using the meteorological observations and global climate data, respectively.

stakes in different periods are shown in Fig. 7. Overall, the model
performs well at stake locations with the RMSE value of 0.003 m
w.e. d1, compared to the mean melt rate of 0.03 m w.e. d1. The
overall correlation coefﬁcient r between measured and modeledablation in different periods is 0.78 (Fig. 7a). In particular, the correlation coefﬁcient r between measured and modeled ablation in
the summer of 2008 is 0.82 (Fig. 7b). However, the calibrated model overestimates the low melt rates beneath the thick debris layer(Fig. 7a). This may be inﬂuenced by the water content of the thickdebris layer. The water content of the debris layer exerts a strong
inﬂuence on the melt rate beneath the debris layer through affecting its thermal properties (Nakawo and Young, 1981). Nakawo andYoung (1981) found that the thermal conductivity of the debris
layer increases with the water content of the layer increasing,resulting in the melt rate beneath the debris layer increasing. They
also found that modeled melt rates for dry thick debris is larger
than those for saturated thick debris. However, considering the effect of the water content must allow the physical properties of thedebris layer to vary with depth, which would require detailed dataon internal debris temperature and water content, and these properties are difﬁcult to determine for an entire glacier. Hence, thisproperty is unlikely to be known in any practical application of
the model presented in previous studies (e.g. Nakawo and Young,
1981; Nicholson and Benn, 2006; Reid and Brock, 2010). In thisstudy we calculate the melt rate to most closely approximate themeasurements without considering the nature of the water content of the debris layer, despite the fact that the debris may be saturated where the layer is thick, and moist where it is thicker, so themeasurements will be expected to plot closer to the modeled results for saturated debris. This pattern suggests that the model systematically overestimates the melt rate in this environment.
Despite this uncertainty, overall, the model can produce reliable
estimates of ablation beneath the debris layer on an annual timescale (Fig. 7) and allows us to consider the effect of debris coveron modifying the catchment mass balance.

Fig. 6. Evaluation of the calibrated model using scatter diagrams of observed versus
modeled monthly runoff (a) and cumulative runoff (b). Qmod, observed data andQmod, global data denote the simulations using the meteorological observations
and global climate data, respectively.

Overall the calibrated model performs well and can simulate
the observations accurately. This indicates that the model has predictive power for mass balance modeling on a catchment scale,assuming that the calibrated model parameter is valid for different

156

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Fig. 8. (a) Calculated speciﬁc mass balance (SMB) showing both the debris-covered
surface and the assumed debris free on glaciers, (b) calculated equilibrium line
altitude (ELA), (c) mean annual temperature (T), and (d) annual precipitation (P).
Light gray shading for the SMB and ELA denotes the standard deviation. Lines in (b),(c), and (d) denote the average values of ELA, T and P over the period 1952–2009.
Fig. 7. Evaluation of the calibrated model using scatter diagrams of observed vsmodeled stake ablation on Hailuogou Glacier in 1982/83, 1990–1994 (a), and thesummer of 2008 (b).

6. Discussion
6.1. Model uncertainty

time periods. In particular, the model is able to reliably estimate
glacier ablation beneath debris, allowing us to consider the inﬂuence of debris on the mass balance of the catchment. We evaluate
the use of global climate data for modeling in Hailuogou catchment
and found that accurate estimation of catchment mass balance ispossible using these data (Fig. 4).
5.3. Mass balance reconstruction
We calculate the long-term mass balance and ELA for Hailuogou
catchment using observed climate data (1989–2009) extended
with global climate data (H08 and Y09) (1951–1988) (Fig. 8). The
ELA is calculated as the altitude where the annual mass balanceis zero. The annual mass balances of Hailuogou catchment are consistently negative (Fig. 8a), with eight positive balance years overthe 58 yrs period (1952–2009), varying between 0.19 and –1.72 m w.e. The mean annual mass balance in the whole catchment
is 0.42 m w.e. for the period 1952–2009. In particular, in recent
years (2001–2009) the mean mass balance is 0.79 m w.e., about
2.3 times the mean value during 1952–2000 (0.35 m w.e.). The
cumulative mass balance of the catchment is 24.3 m w.e. forthe period 1952–2009. In recent years (2001–2009), the cumulative mass balance of the catchment accounts for 30% of the total
cumulative mass balance for the period 1952–2009. The ELA trendof Hailuogou catchment over the period 1952–2009 is shown in
Fig. 8b. The ELA of the glaciers in the catchment varies between
4907 and 5273 m a.s.l. during the study period, with a mean of5068 m a.s.l. The ELA of the glaciers shows a slight negative trendwith a rate of 0.5 m yr1during 1952–2000. In contrast, glaciersshow a signiﬁcant positive trend (12.3 m yr1) of the ELA during
2001–2009 (Fig. 8b). This indicates that the ascending ELA results
in an acceleration of glacier wastage in Hailuogou catchment in recent years.

In highly glacierized catchments, major challenges for mass balance simulation are the lack of available precipitation data and thelarge spatial variability of precipitation due to orographic effects.
Furthermore, the precipitation gradient is difﬁcult to estimate,which makes spatial interpolation challenging. In this study, theprecipitation increase with altitude is set to a ﬁxed gradient thatis optimized by the calibration procedure, minimizing errors between the model results and observations. However, uncertainties
still arise due to the high spatial variability of precipitation, and theprecipitation gradient is optimized to the mean conditions over theentire study period. In particular, when applying such gradient toall years, larger model errors are to be expected for years when
the observations have larger deviations from the mean conditions.
However, the model predictive ability for long-term mass balance
would be substantially reduced without a generalized parameterset for the entire study period. The main aim of this study is to develop a simple approach that considers the signiﬁcance of debris
cover and its inﬂuence for accurate estimation of the long-term
variation and trend in surface mass balance on the catchment
scale. Hence, a generalized parameter set was not unexpected to
be used for the mass balance calculation in Hailuogou catchment.
Glacier area data are generally not available every year. Somestudies assumed the constant glacier area in their mass balancecalculations (e.g. Fischer, 2009; Fujita et al., 2011). Hirabayashi
et al. (2010) presented an approach in which the glacier area perelevation band is updated daily. In our analysis, glacier areas in
1966, 1975, 1994, and 2007 are used. These are correspondingly
assumed to be constant over the periods 1952–1966, 1967–1975,
1976–1994, and 1995–2009. As reported by Liu et al. (2010), theglacier area of the catchment decreased by 2.9 km2 during 1966–
2007 (7.8% of the 1966 total area), which mainly occurred in thelower part of the glaciers. To assess the uncertainties in constantY. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

glacier areas in different periods, the mass balance during the entire analysis period is calculated using the glacier areas in 1966,
1975, 1994, and 2007, respectively. The resulting mass balance
indicates that the speciﬁc mass balances estimated using the glacier areas of 1966, 1975, and 1994 are more negative by 0.06,
0.07, and 0.04 m w.e. (the average of the calculation period),
respectively, compared to those estimated by using the glacierareas in 1975,1994, and 2007. This ﬁnding indicates that changes
in glacier extent during the periods 1952–1966, 1967–1975,
1976–1994, and 1995–2009 have a small impact on the speciﬁcmass balance of Hailuogou catchment.
6.2. Effect of supraglacial debris on glacier mass balance and runoff
Hailuogou catchment is dominated by debris-covered glaciers(Table 1). The total debris-covered area accounts for 39% of the total area of ablation zones of three debris-covered glaciers in thecatchment. Zhang et al. (2011) extensively surveyed the debris
thickness on Hailuogou Glacier, where the debris-covered area accounted for 41% of the total area of the ablation zone. They suggested that the debris thickness varies from just several
millimetres of patchy cover to >1 m, with an inhomogeneous spatial distribution. As suggested by previous studies, dispersed andthin debris enhances ice melt rates through albedo reduction,
whereas debris cover of thickness exceeding a few centimetres reduces ice melt by insulating the ice from atmospheric heat(Østrem, 1959; Nakawo and Young, 1981; Mattson et al., 1993).
Zhang et al. (2011) found that because of the inhomogeneous distribution of debris thickness, about 67% of the ablation area on
Hailuogou Glacier 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 in only 14%of the ablation area; on the other hand, the altitude distributionof the debris thickness further complicated the ablation gradienton the ablation area, making it completely different to that of debris-free glaciers. Therefore, spatial distribution of debris plays acrucial role in the estimation of glacier ablation, and must be considered in the mass balance estimation of Hailuogou catchment. To
assess the effect of the debris cover on the mass balance of Hailuogou catchment during the past few decades, we recalculate themass balance assuming that there is no debris on the three glaciersin the catchment. In accordance with the calculation, we ﬁnd a 12%
increase of glacier melt under the real surface condition during
1952–2009 compared to that with an assumption of no debris cover on the glaciers. As a result, the positive mass balance years account for 35% of all years (58 yrs) without debris cover on theglaciers, but they represent only 14% of all years under the real surface conditions (Fig. 8a). The change in mean mass balance withthe real surface condition during 1952–2009 is twice as much asthat without debris cover on the glaciers (0.21 m w.e.). This indicates that supraglacial debris cover signiﬁcantly accelerates massloss from the glaciers in Hailuogou catchment (Fig. 8a).
Scherler et al. (2011; Fig. 2) reported that widespread debris cover on many Himalayan glaciers reduces their retreat rates, making
their fronts stable, which are therefore unsuitable as indicators of recent climate change. In contrast, the termini of the debris-covered
glaciers in Hailuogou catchment have been retreating at a high rateduring the past few decades (Su et al., 1992; Liu et al., 2010; Zhang
et al., 2010). The terminus of Hailuogou Glacier receded at a meanannual rate of 27.4 m yr1 during 1998–2008 (Zhang et al., 2010),
and Glacier No. 2 receded at a mean annual rate of 30.0 m a1 during
the period 1981–1990 (Su et al., 1992). Based on the glacier outlines
in 1966 and 2007 created by Liu et al. (2010), we calculate the retreat
rates of the termini of seven glaciers in the catchment. The mean annual retreat rate of the termini of three debris-covered glaciers is28.9 m yr1 during 1966–2007, which is signiﬁcantly larger than157

that of four debris-free glaciers (15.6 m yr1). The analysis of thespatial characteristics of the ice melt rates on the ablation area of
Hailuogou Glacier indicated that one of the high melt rates occurs
in the terminus zone (altitudes of 3000–3200 m a.s.l.), an area witha wide distribution of coexisting debris-covered ice, bare ice, icecliffs, and supraglacial ponds (Zhang et al., 2011). Such coexisting
surface that is composed of debris-covered ice, bare ice, ice cliffs
and supraglacial ponds, is also widely distributed in the terminus
zones of other debris-covered glaciers in the catchment (Li and Su,
1996). In these areas, in addition to debris cover and its spatial distribution impacting the ice melt rate, ice cliffs and supraglacial
ponds signiﬁcantly accelerate the ice melt rate (Sakai et al., 2000,
2002). Ice cliffs are generally covered with a thin layer of debris or
dust, meaning that the melt rate is higher than on nearby debriscovered surface (Sakai et al., 2000, 2002). The analysis of a pond heatbudget indicates that the energy absorption at the supraglacial pond
is several times larger than that at the surrounding debris-covered
surface, which results in accelerating the ice melt rate (Sakai et al.,2000). Hence, spatial distribution of the areas with the coexistingdebris-covered ice, bare ice, ice cliffs, and supraglacial ponds leads
to the spatially inhomogeneous ice ablation, which is likely to makea disproportionately large contribution to ablation at the glacier termini. On the other hand, supraglacial ponds can accelerate glacierterminus disintegration by exposing ice faces at the surface (Röhl,
2008) and lowering of the glacier surface (Sakai and Fujita, 2010).
Due to the inhomogeneous ablation, the debris layer becomesincreasingly unstable in these areas (Iwata et al., 2000; Benn et al.,2001), releasing debris via occasional sluicing of debris down steep
slopes. This process redistributes the debris upon the glacier surface,
and then modiﬁes the spatial variability in the ice melt rate, causing
differential surface lowering, development of ice cliffs and evolutionof supraglacial ponds (Iwata et al., 2000; Benn et al., 2001; Sakai
et al., 2000, 2002; Röhl, 2008; Sakai and Fujita, 2010). Furthermore,
the ice velocity is relatively large on debris-covered glaciers in thecatchment (Li and Su, 1996; Zhang et al., 2010). Ice velocities measured on Hailuogou Glacier indicated that the ice velocity at the terminus is 28.8 m yr1 (Zhang et al., 2010), which is signiﬁcantly
larger than that in the Himalayan study of Scherler et al. (2011). Consequently, debris and its spatial distribution in association with ice
cliffs and supraglacial ponds results in the unstable fronts of debris-covered glaciers in this catchment, which accelerates their terminus disintegration and retreat.
The main contributions to the runoff of Hailuogou catchment
are meltwater and rainfall (Liu et al., 2010). Note that glacier meltwater is the main contribution to the increasing runoff of thecatchment during 1994–2002 when the precipitation of the catchment shows a decrease trend (Fig. 8d). As discussed above,
Hailuogou catchment is dominated by debris-covered glaciers,
where the ice melting depends strongly on debris thickness andits spatial continuity. To assess the effect of the debris cover on
the runoff in the catchment, we recalculate the runoff assumingthat there is no debris on three debris-covered glaciers. In accordance with the calculation, we ﬁnd that the mean annual runoff
without debris cover on three glaciers during 1994–2004 decreases
by 6% compared to that with the real surface condition. In addition,
evaporation of rainwater from debris appears to have a small inﬂuence on catchment hydrology and downstream water resources(Brock et al., 2010). It is apparent that the presence of debris coverleads to an increase in catchment runoff through its inﬂuence ofaccelerating surface melt, and its insulation effect is not signiﬁcantin Hailuogou catchment.
6.3. Climate change, mass balance, and ELA
During the past 58 yrs (1952–2009), glaciers in Hailuogou
catchment experience signiﬁcant mass loss (Fig. 8a), especially in158

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Fig. 9. (a and b) Regression relations of summer temperature vs the speciﬁc mass balance (SMB) and the ELA, and (c and d) regression relations of summer precipitation
versus the SMB and the ELA for the period 1952–2009. Here, R2 is the coefﬁcient of determination.

Table 3Mass balance (MB; m w.e.), summer temperature (Temp.; °C), and precipitation(Prec.; mm), and their contributions to the glacier mass balance in each decade
estimated based on the regressions shown in Fig. 9.
Temp.

1950s
1960s
1970s
1980s
1990s
2000s

9.09
9.19
9.02
9.05
9.46
9.82

Prec.

1473.11450.91449.71513.91541.21508.9

MB

0.41
0.51
0.36
0.20
0.21
0.78

Variable contribution (%)
Temp.

Prec.

19.6
35.9
15.2
–
135
97.5

80.4
64.1
115.2
–
35
2.5

the last 10 yrs. These ﬁndings agree with those of Liu et al. (2010)
and Zhang et al. (2010), who reported an accelerated shrinkage ofglaciers in the last few decades based on topographic maps, remotesensing data, and the GAEORS meteorological data. The glaciers inthis catchment are the summer-accumulation type, which is moresensitive to climate change (Shi and Liu, 2000). As shown in Fig. 8c,
the air temperature of the catchment shows a signiﬁcant increasing trend with a rate of 0.13 °C (10 yr)1 over the period 1952–2009, especially in the last 20 years when the rate was 0.27 °C(10 yr)1. In addition, the air temperature at the beginning of themelt season (April) and at its end (October) rose by 1.6 and
0.9 °C, respectively, for the period 1998–2004 compared to thosefor the period 1988–1997 (Zhang et al., 2010). This possibly prolongs the melt period and causes a signiﬁcant decrease in snow

accumulation, and an additional decrease in surface albedo due
to a diminished summer snowfall that prevents the absorption ofsolar radiation and snow melt during the melt season (Fujita,2008). On the other hand, the annual precipitation shows a slightincreasing trend during 1952–2009 (Fig. 8d). As shown in Fig. 8d,
annual precipitation is relatively low in the 1950s, 1960s, and1970s, but high in the 1980s and 1990s. After 2000 annual precipitation decreases compared to that in the 1990s. In addition, theprecipitation from July to September decreased by 24% for the period 1998–2004 compared to that for the period 1988–97, especially in July and September (Zhang et al., 2000). This change
possibly has a signiﬁcant effect on glacier accumulation in thecatchment.
To quantify the effect of climate change on glaciers in the catchment, we examine the relationships among summer (May–
October) temperature, precipitation, mass balance and ELA. Figs.
9a and b shows the relations between summer temperature andthe mass balance and ELA, respectively. Summer temperature isnegatively correlated with the mass balance and positively associated with the ELA. The regression results in Fig. 9a and b suggestthat the 1 °C summer temperature change results in 0.644 m w.e.
glacier mass loss, and 135.3 m ELA ascending, respectively. Similarly, summer precipitation is positively correlated with the massbalance and negatively associated with the ELA (Fig. 9c and d).
The regression results in Fig. 9c and d suggest that the 100 mmsummer precipitation change leads to 0.25 m w.e. mass accumulation, and 36.7 m ELA descending, respectively. Compared to theprecipitation-mass balance/ELA relations, those between air

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

temperature and mass balance/ELA are stronger, especially in recent years. Therefore, we can investigate which variables have controlled the glacier mass balance in each decade using the 1980s as
a reference decade. As shown in Table 3, summer precipitation issigniﬁcantly reduced in the 1950s, 1960s, and 1970s compared toin the 1980s, and summer temperature shows a slight variation
over the same period. Based on the regression results mentioned
above, a marked reduction in summer precipitation that was observed before the 1980s appears to be the main cause of the negative mass balance. After the 1980s, the summer temperature showsa signiﬁcant increasing trend, and the precipitation shows a slightincreasing trend in the 1990s and a slight decreasing trend in the
2000s (Table 3). Increasing summer temperature is the main causeof rapid wasting of the glacier mass in Hailuogou catchment (Table3). In particular, the contribution of the rise in temperature to glacier mass loss increases to 97.5% in the 2000s (Table 3). This indicates that temperature change controls variation in glacier massbalance in Hailuogou catchment over the last 20 yrs.

7. Conclusions
A surface energy–mass balance model is implemented on acatchment scale, and the signiﬁcance of debris and its inﬂuenceon glacier melt rates is taken into account. The model is used toreconstruct the mass balance time series of glaciers in Hailuogou
catchment in the south-eastern Tibetan Plateau, where there arethree debris-covered and four debris-free monsoonal maritime
glaciers. Global climate data (1952–1988) and observations(1989–2009) are used to drive the model. The global climate data
are validated against long-term meteorological observations from
Hailuogou catchment and correspond sufﬁciently well with surface
observations to be used for mass balance modeling. The model iscalibrated against runoff observations in the catchment using multiple assessment criteria to minimize the error between simulations and observations. It is veriﬁed by comparing the model
results to the observed runoff and the ablation measured at stakes
in different periods. Overall, the calibrated model performs well,which indicates that the model has predictive power for mass balance modeling on a catchment scale, assuming that the calibrated
model parameter is valid for different time periods. In particular,
the model is able to reliably estimate glacier ablation beneath debris, allowing us to consider the effect of debris on the mass balanceof the catchment.
From the reconstructed time series, a mass loss of 24.3 m w.e.
is estimated in Hailuogou catchment over the period 1952–2009.
The annual mass balances are consistently negative over the
58-yr period with eight positive balance years. In particular, glaciers in the catchment show markedly rapid wastage with an associated increasing trend of the ELA in recent years (2001–2009).
Regression analyses show that the relations between summer temperature and mass balance/ELA are much stronger relative to theprecipitation-mass balance/ELA relations. On the basis of these
regression relations, we found that a marked reduction in summer
precipitation is the main cause of the negative mass balance before
the 1980s, whereas signiﬁcant warming in the summer is the maincause of rapid wasting of the glacier mass in the catchment in recent decades. In addition, recalculation with an assumption of no
debris cover on the glaciers in the catchment indicates that thepresence of debris results in the 12% increase of glacier melt during
1952–2009, and then signiﬁcantly accelerates glacier mass loss.
Furthermore, debris and its spatial distribution in association withice cliffs and supraglacial ponds results in the unstable fronts ofdebris-covered glaciers in Hailuogou catchment.
Debris-covered glaciers are not restricted to the Tibetan Plateau,
but are common in many other mountain ranges around the world.

159

For modeling mass balance, predicting runoff and response to
climate change, debris cover and its inﬂuence on glacier melt rates
should be added to analysis that has determined glacier mass balances. Using the energy–mass balance model presented in thisstudy, we can accurately estimate the long-term variation andtrend in surface mass balance on a catchment scale and examine
the signiﬁcance of debris cover and its inﬂuence on glacier melts.
The model can also output catchment runoff. The global climate
data used in this study can be used to drive such a model. Thismeans that the model can be applied to other catchments in the
Tibetan Plateau, where debris-covered glaciers exist, to investigate
glacier response to climate change and the impact of global warming on the regional water cycle, even in regions that are not routinely monitored because of the harsh climatic conditions andremoteness.

Acknowledgements
We thank to Trend Biggs, Koji Fujita and the anonymous reviewer for their extensive and very helpful comments and suggestions. We also thank the Gongga Alpine Ecosystem Observation
and Research Station of the Chinese Ecological Research Network
for providing the meteorological data. This research was supportedby the Global Change Research Program of China (2010CB951401),
the National Natural Science Foundation of China (40701032,
40801030), and the Funding Program for next generation worldleading researchers, JSPS and CREST project of Japan Science andTechnology Agency (JST).

References
Aizen, V.B., Aizen, E.M., 1994. Regime and mass-energy exchange of subtropicallatitude glaciers under monsoon climatic conditions: Gongga Shan, Sichuan.
China. J. Mt. Res. Dev. 14 (2), 101–118.
Arnold, N.S., Willis, C., Sharp, M.J., Richards, K.S., Lawson, W.J., 1996. A distributed
surface energy balance model for a small valley glacier: I. Development andtesting for Haut Glacier d’Arolla, Valais, Switzerland. J. Glaciol. 42 (140), 77–89.
Barry, R.G., 2006. The status of research on glaciers and global glacier recession: areview. Prog. Phys. Geogr. 30, 285–306. http://dx.doi.org/10.1191/0309133306
pp 0309133478ra.
Benn, D.I., Wiseman, S., Hands, K.A., 2001. Growth and drainage of supraglacial lakes
on the debris-mantled Ngozumpa Glacier, Khumbu Himal. Nepal. J. Glaciol. 47(159), 626–638.
Bougamont, M., Bamber, J.L., Greuell, W., 2005. A surface mass balance model for
the Greenland ice sheet. J. Geophys. Res. 110, F04018. http://dx.doi.org/
10.1029/2005JF000348.
Braithwaite, R.J., Zhang, Y., 1999. Modelling changes in glacier mass balance that
may occur as a result of climate changes. Geogr. Ann. 81A, 489–496.
Braun, L.N., Grabs, W., Rana, B., 1993. Application of conceptual precipitation runoff
model in the Langtang Khola Basin, Nepal Himalaya. In: Young, G.J. (Eds.), Snowand Glacier Hydrology. IAHS Publication, Proceedings of the Kathmandu
Symposium, November 1992, pp. 221–237.
Brock, B.W., Mihalcea, C., Kirkbride, M.P., Diolaiuti, G., Cutler, M.E.J., Smiraglia, C.,
2010. Meteorology and surface energy ﬂuxes in the 2005–2007 ablation seasonsat the Miage debris-covered glacier, Mont Blanc Massif, Italian Alps. J. Geophys.
Res. 115, D09106. http://dx.doi.org/10.1029/2009JD013224.
Cao, Z., Cheng, G., 1994. Preliminary analyses of hydrological characteristics ofHailuogou Glacier on the eastern slope of the Gongga Mountain. In: Xie, Z.,
Kotlyakov, V.M. (Eds.), Glaciers and Environment in the Qinhai-xizang (Tibet)
Plateau (I) – The Gongga Mountain. Science Press, Beijing and New York, pp.
143–156.
Cheng, G., 1996. Exploration of precipitation features on extra-high zone of Mt
Gongga. Mt. Res. 14 (3), 177–182 (In Chinese with English abstract).
Dyurgerov, M.B., Meier, M.F., 2005. Glaciers and the changing earth system: a 2004
snapshot. Occasional Paper 58, Institute of Arctic and Alpine Research,
University of Colorado at Boulder, Boulder, Colorado, USA.
Fan, Y., van den Dool, H.V.D., 2008. A global monthly land surface air temperature
analysis for 1948-present. J. Geophys. Res. 113, D01103. http://dx.doi.org/
10.1029/2007JD008470.
Fischer, A., 2009. Glaciers and climate change: interpretation of 50 years of direct
mass balance of Hintereisferner. Global Planet. Change 71, 13–26. http://dx.doi.org/10.1016/j.gloplacha.2009.11.014.
Fujita, K., 2007. Effect of dust event timing on glacier runoff: sensitivity analysis fora Tibetan glacier. Hydrol. Process. 21 (21), 2892–2896. http://dx.doi.org/
10.1002/hyp. 6504.

160

Y. Zhang et al. / Journal of Hydrology 444–445 (2012) 146–160

Fujita, K., 2008. Effect of precipitation seasonality on climatic sensitivity of glaciermass balance. Earth Planet. Sci. Lett. 276, 14–19. http://dx.doi.org/10.1016/
j.epsl.2008.08.028.
Fujita, K., Ageta, Y., 2000. Effect of summer accumulation on glacier mass balance on
the Tibetan Plateau revealed by mass–balance model. J. Glaciol. 46 (153), 244–
252.
Fujita, K., Nuimura, T., 2011. Spatially heterogeneous wastage of Himalayan
glaciers.
PNAS
108
(34),
14011–14014.
http://dx.doi.org/10.1073/
pnas.1106242108.
Fujita, K., Seko, K., Ageta, Y., Pu, J., Yao, T., 1996. Superimposed ice in glacier massbalance on the Tibetan Plateau. J. Glaciol. 42 (142), 454–460.
Fujita, K., Takeuchi, N., Nikitin, S.A., Surazakov, A.B., Okamoto, S., Aizen, V.B., Kubota,J., 2011. Favorable climatic regime for maintaining the present-day geometry of
the Gregoriev Glacier, Inner Tien Shan. TC 5, 539–549. http://dx.doi.org/
10.5194/tc-5-539-2011.
Fukushima, Y., Kawashima, K., Suzuki, M., Ohta, T., Motoyama, H., Kubota, H., Yamada,T., Bajracharya, O.R., 1987. Runoff characteristics in three glacier-covered
watersheds of Langtang Valley, Nepal Himalaya. Bull. Glaciol. Res. 5, 11–18.
Hirabayashi, Y., Kanae, S., Emori, S., Oki, T., Kimoto, M., 2008. Global projections ofchanging risks of ﬂoods and droughts in a changing climate. Hydrolog. Sci. J. 53(4), 754–772.
Hirabayashi, Y., Döll, P., Kanae, S., 2010. Global-scale modeling of glacier massbalances for water resources assessments: Glacier mass changes between 1948
and
2006.
J.
Hydrol.
390,
245–256.
http://dx.doi.org/10.1016/
j.jhydrol.2010.07.001.
Hock, R., 1999. A distributed temperature-index ice- and snowmelt model including
potential direct solar radiation. J. Glaciol. 45, 101–111.
Hock, R., Holmgren, B., 2005. A distributed surface energy-balance model forcomplex topography and its application to Storglaciären, Sweden. J. Glaciol. 51(172), 25–36.
Hock, R., Radić, V., Woul, M.D., 2007. Climate sensitivity of Storglaciären, Sweden:
an intercomparison of mass-balance models using ERA-40 re-analysis andregional climate model data. Ann. Glaciol. 46, 342–348.
Immerzeel, W.W., van Beek, L.P.H., Konz, M., Shrestha, A.B., Bierkens, M.F.P., 2012.
Hydrological response to climate change in a glacierized catchment in the
Himalayas. Clim. Change 110, 721–736. http://dx.doi.org/10.1007/s10584-0110143-4.
Iwata, S., Aoki, T., Kadota, T., Seko, K., Yamaguchi, S., 2000. Morphological evolutionof the debris cover on Khumbu Glacier, Nepal, between 1978 and 1995. In:
Raymond, C.F., Nakawo, M., Fountain, A. (Eds.), Debris-Covered Glaciers. IAHS
Publication, Wallingford, UK, Proceedings from a Workshop Held at Seattle, WA,
USA, pp. 3–11.
Jarvis, A., Reuter, H.I., Nelson, A., Guevara, E., 2008. Hole-ﬁlled seamless SRTM data
V4, International Centre for Tropical Agriculture (CIAT) <http://
srtm.csi.cgiar.org/> (accessed: 12.03.11).
Kang, E., Cheng, G., Lan, Y., Jin, H., 1999. A model for simulating the response of
runoff from the mountainous watersheds of inland river basins in the arid zoneof northwest China to climatic changes. Sci. China Ser. D 42, 52–63.
Klok, E.J., Oerlemans, J., 2002. Model study of the spatial distribution of the energy
and mass balance of Morteratschgletscher, Switzerland. J. Glaciol. 48 (163),
505–518.
Kondo, J. (Ed.), 1994. Meteorology of Water Environment. Asakura, Tokyo (InJapanese).
Lemke, P., Ren, J., Alley, R.B., Allison, I., Carrasco, J., Flato, G., Fujii, Y., Kaser, G., Mote,P., Thomas, R.H., Zhang, T., 2007. Observations: changes in snow, ice and frozen
ground. In: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K.B.,
Tignor, M., Miller, H.L. (Eds.), Climate Change 2007: The Physical Science Basis.
Cambridge University Press, Cambridge, United Kingdom and New York, NY,
USA, p. 356.
Li, J., Su, Z., 1996. Glaciers in the Hengduan Mountains. Science Press, Beijing (InChinese, pp. 1–110).
Liu, S., Xie, Z., Su, Z., Song, G., Wang, N., 1994. Accumulation, ablation and ice
formation of Hailuogou Glacier in the Gongga Mountain. In: Xie, Z.a.V.M.K. (Ed.),
Glaciers and environment in the Qinhai-xizang (Tibet) Plateau (I) –The Gongga
Mountain. Science Press, Beijing and New York, pp. 65–80.
Liu, S., Zhang, Y., Zhang, Y., Ding, Y., 2009. Estimation of glacier runoff and future
trends in the Yangtze River source region. China. J. Glaciol. 55 (190), 353–362.
Liu, Q., Liu, S., Zhang, Y., Wang, X., Zhang, Y., Guo, W., Xu, J., 2010. Recent shrinkage
and hydrological response of Hailuogou glacier, a monsoon temperate glacieron the east slope of Mount Gongga. China. J. Glaciol. 56 (196), 215–224.
Matsuda, Y., Fujita, K., Ageta, Y., Sakai, A., 2006. Estimation of atmospherictransmissivity of solar radiation from precipitation in the Himalayas andTibetan Plateau. Ann. Glaciol. 43, 344–350.
Mattson, L.E., Gardner, J.S., Young, G.J., 1993. Ablation on debris covered glaciers: an
example from the Rakhiot Glacier, Punjab, Himalaya. In: Young, G.J. (Eds.), Snowand Glacier Hydrology. IAHS Publication, Proceedings of the Kathmandu
Symposium, November 1992, pp. 289–296.

Mihalcea, C. et al., 2006. Ice ablation and meteorological conditions on the debriscovered area of Baltoro glacier, Karakoram, Pakistan. Ann. Glaciol. 43,
292–300.
Mitchell, T.D., Jones, P.D., 2005. An improved method of constructing a database ofmonthly climate observations and associated high-resolution grids. Int. J.
Climatol. 25, 693–712. http://dx.doi.org/10.1002/joc.1181.
Motoyama, H., 1990. Simulation of seasonal snowcover based on air temperatureand precipitation. J. Appl. Meteorol. 29 (11), 1104–1110.
Nakawo, M., Young, G.J., 1981. Field experiments to determine the effect of a debris
layer on ablation of glacier ice. Ann. Glaciol. 2, 85–91.
Nakawo, M., Young, G.J., 1982. Estimate of glacier ablation under a debris layer fromsurface temperature and meteorological variables. J. Glaciol. 28 (98), 29–34.
Nash, J.E., Sutcliffe, J.V., 1970. River ﬂow forecasting through conceptual models.
Part 1. A discussion of principles. J. Hydrol. 10 (3), 282–290.
Nicholson, L., Benn, D.I., 2006. Calculating ice melt beneath a debris layer using
meteorological data. J. Glaciol. 52 (178), 463–470.
Oerlemans, J., Fortuin, J.P.F., 1992. Sensitivity of glaciers and small ice caps togreenhousewarming.
Science258,
115–117,
doi:110.1126/
science.1258.5079.1115.
Ohmura, A., 2001. Physical basis for the temperature-based melt-index method. J.
Appl. Meteorol. 40, 753–761.
Ohmura, A., 2006. Changes in mountain glaciers and ice caps during the 20th
century. Ann. Glaciol. 43, 361–368.
Østrem, G., 1959. Ice melting under a thin layer of moraine and the existence of ice
cores in moraine ridges. Geogr. Ann. 41, 228–230.
Pu, J. (Ed.), 1994. Glacier Inventory of China VIII (The Changjiang River Drainage
Basin). Gansu Culture Press, Lanzhou, pp. 117–129 (In Chinese).
Reichert, B.K., Bengtsson, L., Oerlemans, J., 2001. Midlatitude forcing mechanisms
for glacier mass balance investigated using general circulation models. J.
Climatol.
14,
3767–3784.
http://dx.doi.org/10.1175/1520-0442(2001)
3014<3767:MFMFGM>3762.3760.CO;3762.
Reid, T.D., Brock, B.W., 2010. An energy-balance model for debris-covered glaciers
including heat conduction through the debris layer. J. Glaciol. 56 (199), 903–
916.
Rignot, E., Rivera, A., Casassa, G., 2003. Contribution of the Patagonia Iceﬁelds ofSouth America to sea level rise. Science 302, 434–437. http://dx.doi.org/
10.1126/science.1087393.
Röhl, K., 2008. Characteristics and evolution of supraglacial ponds on debriscovered Tasman Glacier. New Zealand. J. Glaciol. 54 (188), 867–880.
Sakai, A., Fujita, K., 2010. Formation conditions of supraglacial lakes on debriscovered glaciers in the Himalaya. J. Glaciol. 56 (195), 177–181.
Sakai, A., Takeuchi, N., Fujita, K., Nakawo, M., 2000. Role of supraglacial ponds in theablation processes of a debris-covered glacier in the Nepal Himalayas. In:
Raymond, C.F., Nakawo, M., Fountain, A. (Eds.), Debris-Covered Glaciers. IAHS
Publication, Wallingford, UK, Proceedings from a Workshop held at Seattle, WA,
USA, pp. 119–130.
Sakai, A., Nakawo, M., Fujita, K., 2002. Distribution characteristics and energybalance of ice cliffs on debris-covered glaciers, Nepal Himalaya. AAAR 34 (1),
12–19.
Sakai, A., Fujita, K., Nakawo, M., Yao, T., 2009. Simpliﬁcation of heat balancecalculation and its application to the glacier runoff from the July 1st Glacier innorthwest China since the 1930s. Hydrol. Process. 23 (4), 585–596. http://dx.doi.org/10.1002/hyp. 7187.
Scherler, D., Bookhagen, B., Strecker, M.R., 2011. Spatially variable response of
Himalayan glaciers to climate change affected by debris cover. Nat. Geosci. 4,
156–159. http://dx.doi.org/10.1038/ngeo1068.
Shi, Y., Liu, S., 2000. Estimation on the response of glaciers in China to the global
warming in the 21st century. Chinese Sci. Bull. 45, 668–672.
Su, Z., Liu, S., Wang, N., Shi, A., 1992. Recent ﬂuctuations of glaciers in the Gongga
Mountains. Ann. Glaciol. 16, 163–167.
Su, Z., Song, G., Cao, Z., 1996. Maritime characteristics of Hailuogou Glacier in the
Gongga Mountains. J. Glaciol. Geocryol. 18, 51–59 (In Chinese with English
abstract).
Yasunari, T.J., Koster, R.D., Lau, K.-M., Aoki, T., Sud, Y.C., Yamazaki, T., Motoyoshi, H.,
Kodama, Y., 2011. Inﬂuence of dust and black carbon on the snow albedo in theNASA Goddard Earth Observing System version 5 land surface model. J.
Geophys. Res. 116, F03017. http://dx.doi.org/10.1029/2010JF001939.
Yatagai, A. et al., 2009. A 44-year daily gridded precipitation dataset for Asia basedon a dense network of rain gauges. SOLA 5, 137–140. http://dx.doi.org/10.2151/
sola.2009-035.
Zhang, Y., Fujita, K., Liu, S., Liu, Q., Wang, X., 2010. Multi-decadal ice velocity andelevation changes of a monsoonal maritime glacier, Hailuogou Glacier. China. J.
Glaciol. 56 (195), 65–74.
Zhang, Y., Fujita, K., Liu, S., Liu, Q., Nuimura, T., 2011. Distribution of debris
thickness and its effect on ice melt at Hailuogou Glacier, southeastern
Tibetan Plateau, using in situ surveys and ASTER imagery. J. Glaciol. 57 (206),
1147–1157.