Hydrological Research Letters 7(1), 6–11 (2013)
Published online in J-STAGE (www.jstage.jst.go.jp/browse/HRL). doi: 10.3178/HRL.7.6

Supplement Information

Supplement Information S1. Model sensitivity to precipitation lapse rate
In HYOGA2, precipitation is assumed to be constant for all elevation bands ineach grid cell. Nevertheless, in order to examine the effect of precipitation lapse rate,
a simple sensitivity test assuming a fixed precipitation lapse rate (a 10% increase per
100 m, similar to the value used in Radić and Hock (2011)) was applied to glaciers
where observed mass balance can be obtained from DM05. Supplement Figure S1(d)shows the result of this sensitivity test: correlation coefficients between the computed
mass balance with the fixed precipitation lapse rate and the observed mass balance of
the available stations in DM05. Supplement Figure S1(d) is similar to Supplement
Figure S1(c) which is the result of “without precipitation lapse rate”; for example,
these two figures are similar in terms of the mean and standard deviation ofcorrelation coefficients. As such, the result indicated that the modeled mass balance is
not sensitive to this parameter in the retrospective simulation.
The same sensitivity test was then extended to a future warmer climate derived
from a single GCM simulation, MIROC5. Thirty four stations among the 114 stationsin Supplement Figure S1, where the model was well calibrated and glacier remainedin 2100, were selected for the test. On average, the difference in cumulative massbalance since 1948 between the control simulation and the simulation with the fixed
precipitation lapse rate was 4% in current climate (~2005) and 10% at the end of the
future climate (~2100). Nonetheless, as was described in the main text, we basicallydid not use precipitation lapse rates in the simulations to avoid the incorporation of anadditional uncertain parameter. In reality, precipitation lapse rates on a global scale
are unknown. We would argue that the incorporation of adequate precipitation lapse

1rate is a subject of future research and is beyond the scope of this letter. Further
studies are expected to examine the effect of precipitation lapse rate in a future
warming climate and to set it adequately on a global scale. A future study is also
expected to investigate the physical reason behind the effect of precipitation lapse rateon mass balance estimation.

Supplement Information S2. Model Calibration and Initialization
The calibration parameters of HYOGA2 are the degree-day factors (DDFs) for
ice and snow. The globally distributed DDFs for HYOGA2 were obtained following amethod proposed by Hirabayashi et al. (2010). The individual glaciers obtained fromthe RGI inventory were first aggregated into one large glacier in 0.5° grid cells. Then,the DDFs for each cell were calibrated until yielding the maximum agreement withthe cell-specific long-term (1948–1980) average of total glacier mass balance
observed at 295 measurement sites. The mass balance observations, glacier areas andlocations were available from Dyurgerov and Meier (2005) (DM05). If there were no
mass balance observations within five grid cells from the target cell, and the glaciergrid cell was located within one of the five regions with available data, the averageregional mass balances of Dyurgerov and Meier (1997) and Serreze et al. (2000) wereused for calibration. These regional balances were obtained as the area-weighted
means of the observed mass balances. In total, 2545 out of the 3120 aggregatedglaciers at 0.5° grid cells were calibrated against these cell-specific or regional
estimates of long-term glacier mass balance. For the rest of the 575 grid cells in SouthAmerica, Australasia, Africa, the sub-Antarctic islands, and the small Arctic islands,
the long-term average mass balance in 1948–1980 was assumed to be zero. For
estimating the DDFs, the upper and lower limits were set at 1 to 4 (mm °C-1 day-1) for2

snow and 4 and 20 (mm °C-1 day-1) for ice. These ranges were subjectively decided
based on published papers with the objective of obtaining reasonable mass balancechanges compared to the observed values within this model structure. In someglaciers, therefore, these DDF ranges may not reflect a similar range of locallyobserved DDFs. Because the calibrated DDFs, as the only calibration parameter,
might have compensated for all potential errors in the model simulations, model
improvements such as consideration of debris effects could improve the calibration
results. Inclusion of such additional glacier processes along with collection of locallyderived DDFs may improve the calibrated DDF values and initial ranges.
In the original HYOGA, the initial glacier area was considered to be a calibration
parameter. It was calibrated until the calculated glacier area was close to the observed
glacier area at the time of observation. Because HYOGA2 estimated the mass change
in each individual glacier, the calibration of the initial glacier area led to unstable anderroneous simulations in relatively small glaciers. Assuming that the change in glacierarea from the beginning of the model simulation (1948) to the observation year was
negligible for the future area change of all glaciers, the glacier area recorded in the
RGI was taken as an initial area on January 1, 1948. For glaciers located at relatively
low altitudes, where a warming climate has been recorded in the past few decades, theretreating trends of glaciers may already have been obvious in the 20th century. For
these glaciers, it was assumed that the change in area and associated change involume in the past were smaller than the potential change in the future climate. Hence,
the area feedback to altitude was not calculated for the past. For initialization of thesnow amount, the state of snow thickness in each elevation band obtained at the endof the 33-year (1948–1980) simulation was used as an initial condition for the snow
amount in the simulation used for analysis.

3Supplement Information S3. Glacier Inventory
The glacier inventory used in this study was the version 2.0 of the Randolph
Glacier Inventory (RGI, http://www.glims.org/RGI/randolph.html), which wasprocessed by Marzeion et al. (2012). The location (longitude and latitude), maximum
and minimum altitude, and glacier area could be obtained from the inventory. Only
glaciers with complete data on location and area (257,325, with a total area of
515,706 km2) were selected from the RGI (Figure 1). The number of selected MGIlarger than 2 km2 where climate input was available was 21,304, with a total area of
463,410 km2. These glaciers covered ~90% of the total glacier area in the inventory.
In addition, 2,930 LT2s (MGI < 2 km2 in area) with a total area of 898 km2, wereselected from the RGI to represent all LT2s at each 0.5° grid and were considered forcomputation.
For glaciers with the median altitude available in the inventory but without
maximum or minimum altitudes, the maximum or minimum altitude was calculated
by assuming that the median altitude was one-third higher than the minimum altitude,
as suggested by Hirabayashi et al. (2010). For glaciers without a median altitude, butwith either a maximum or minimum altitude, the other value was calculated byassuming that the difference between the maximum and minimum altitudes was 1500m (determined subjectively from glaciers with available data). In addition, for glacierswithout data for both maximum and minimum altitudes, the maximum sample valueof the 30 arc second (approximately 900 m) elevation data of the Shuttle Radar
Topography Mission (SRTM) within each 0.5° grid was assumed to be the maximum
altitude of the glacier. The minimum altitude was taken to be 1500 m lower than the
maximum altitude. Finally, the maximum and minimum altitudes for each glacier4

were adjusted such that they did not exceed the range of maximum and minimum
altitudes of SRTM data within the corresponding 0.5° grid.
For the glaciers in the inventories without data for maximum glacier length Lmax,
the Lmax (m) was calculated from the maximum Zmax (m) and minimum Zmin (m)elevation as
Lmax = K (Zmax – Zmin).

(S1)

The constant ratio K (0.004188) was calculated from glaciers where all three data
were available in the extended World Glacier Inventory (WGI-XF; Cogley, 2010). To
avoid potential model instability, the minimum value of Lmax was set to 25 m, i.e., theminimum length that could be resolved in the 50 m vertical model grid.
As the RGI does not provide the date of observation for individual glaciers, theinitial glacier area in the HYOGA2 historical simulation was assumed to be the sameas the RGI glacier area. Because most of the surveys or images used to define theoutline of RGI were dated after 2000, it was assumed that the impact of the initial
glacier area on future mass balance simulations was minimal.
The mass balance estimation of MGI in HYOGA2 was sensitive to the input
climate forcing, especially to surface temperature. This study used data from a global
database, which may include some grids with unrealistic surface temperatures,
resulting in erroneous mass balance estimation. This is apparent in High Mountain
Asia, where the gauge observations for the correction of global climate datasets islimited. To identify such grids and exclude them from the analysis, upper and lower
limits of annual mass balance were set at 2 and -5 m in the retrospective simulations,
respectively. This range was subjectively decided from the available observed masschanges. We realized that real mass changes in glaciers may be inside or outside thisrange. However, for the purpose of global modeling, we applied this threshold to5exclude potential erroneous calculated mass changes caused by the model. Glaciersshowing erroneous mass balance calculations (before 1960) outside the defined limits
were excluded in the analysis. In total, 2,104 glaciers from the RGI were excluded
from the regional mass balance analysis. The total area of the excluded glaciers was
66,472 km2. We assumed that the exclusion of these glaciers would not significantly
alter the regionally averaged total mass estimates in the model. The ratio of the initial
glacier areas calculated in the model (in total 450,505 km2) and the areas of all thesmall glaciers and ice caps registered in the RGI (in total 515,706 km2) at each 0.5°
grid cell was used to scale up the results to estimate the global volume loss fromglaciers.

Supplement Information S4 . Extension of Observation-based Climate Data
The retrospective simulations for the 64-year period from 1948 to 2011 inclusiveof HYOGA2 were forced by an extended version of the observation-based global 0.5°
gridded dataset of daily precipitation and near-surface temperature (H08; Hirabayashi
et al., 2005, 2008a). The daily temperature of H08 was based on monthly temperatureand monthly diurnal temperature ranges of the Climate Research Unit version TS 2.1
data for 1948–2002 (CRU; Mitchell and Jones, 2005). For extended periods from2003 to 2011, the monthly temperature of Fan and van den Dool (2008) was scaled
using the mean ratio of monthly temperatures from the CRU and Fan and van den
Dool (2008), from 1986 to 2002. The monthly diurnal temperature range of the CRU
TS3.1 (2008–2009) and monthly means of maximum and minimum temperatures of
the Global Telecommunications System (GTS, 2010–2011) were used for extension of
the monthly diurnal temperature range. The daily precipitation data of H08 was basedon version 6 of the monthly precipitation product from the Global Precipitation

6

Climatology Centre (Rudolf and Rubel, 2005 for 1948–2010), and the daily
precipitation product from the GTS (for obtaining the monthly mean for 2011 andwet-day number after 2003) and APHRODITE V1003R1 (Yatagai et al., 2009).
Hirabayashi et al. (2005) described the disaggregation of monthly climatic variables
into a daily time series using a stochastic weather generator, and Hirabayashi et al.
(2008b) explained the under-catch correction for precipitation data. Daily observed
values of APHRODITE with correction of the gauge under-catch bias were used
where and when they were available.

Supplement Information S5. Bias correction of the GCM Data
Because of unavoidable overall biases between GCM simulations and observed
data (e.g., temperature, precipitation), GCM data should be corrected before being
used as inputs in assessment models such as glacier models. It is necessary to bring
the statistical properties of present-day simulations in line with observations and touse this information to correct future projections (Seneviratne et al., 2012). The
original GCM data were first disaggregated onto a 0.5° grid by linear interpolationfrom four neighboring grid cells. Then, the biases in the daily time series of futureland-surface temperature and precipitation were corrected.
The bias correction method explained by Watanabe et al. (2012) was used tocorrect the GCM-simulated temperature and precipitation. The GCM simulation data
(30 years from 1950 to 1980 at 0.5° spatial resolution) were compared against
observation-based daily climate data (H08; Hirabayashi et al., 2008a) to estimate thebiases. The corrections were then conducted in two steps following the method ofHaerter et al. (2011), in which the monthly and daily variances were corrected
simultaneously. First, the monthly bias was corrected by comparison of the monthly

7

GCM and observation data. The daily variability was then corrected by calculating thedifference between the daily time series and monthly means of the GCM data andobservations.
During the bias correction, statistical parameters related to monthly temperatureand precipitation, instead of monthly values, were adjusted. A brief introduction of thebias correction methods for temperature and precipitation are presented here.
The monthly bias correction method for temperature can be expressed as
_

_

_

_

xmon,cor = xmon,obs + xmon,21c – xmon,20c(S2)σmon,cor = σmon,21c σmon,obs (σmon,20c)-1

(S3)

_

_

xmon,cor = f(P(xmon,21c; xmon,21c, σmon,21c); xmon,cor, σmon,cor)

(S4)

_

where xmon,cor and σmon,cor are the corrected mean and standard deviation of monthly__

means of GCM temperature, respectively, xmon,21c and xmon,20c are the means of future_

and current GCM temperature, respectively, xmon,obs is the mean of the observationbased gridded temperature product from H08, σmon,21c and σmon,20c are the standard
deviation for the future and current GCM temperature, and σmon,obs is the standard
deviation of H08 temperature. P(x;a,b) represents a cumulative distribution function(CDF) of the data x, which are normally distributed with mean a and standard
deviation b, and f(x;a,b) represents an inverse function of the CDF P(x;a,b).
The daily time series was then obtained from the bias-corrected monthly time series
as_xday,cor = σday,cor (σday,21c)-1 (xday,21c – xday,21c) + xmon,cor(S5)

where σday,21c is the standard deviation of the daily time series of future GCM
temperature, σday,cor is the standard deviation of the daily time series of corrected8temperature obtained by applying the equidistant CDF matching method (EDCDF
method; Li et al., 2010) to the standard deviations of daily time series ofcurrent/future GCM and observed temperatures. The equidistant CDF matching
method was originally applied to monthly climate data.
The bias correction for monthly precipitation was performed by fitting the
precipitation statistics to a two-parameter gamma distribution, expressed as
_

_

_

_

_

_

kmon,cor = kmon,21c kmon,obs (kmon,20c)-1θmon,cor = θmon,21c θmon,obs (θmon,20c)-1

(S6)

(S7)

where k and θ are the scaling and shape parameters of the gamma distribution,
respectively. The subscript notations are the same as those used for temperature(Equations (S2) through (S4)). Once the bias in monthly GCM precipitation wascorrected, the EDCDF method, as in the case of temperature, was used to compute thebias-corrected daily time series of the precipitation.

Supplement Information S6. Model validation
The old and new versions of the glacier model HYOGA were validated using
observation-based mass balance of glaciers. Figure S1 compares the correlation
coefficient of the computed mass balance and available observations. The mean andstandard deviation of the correlation coefficients among the available stations of
DM05 were 0.47±0.24 for HYOGA, 0.48±0.25 for HYOGA2 with calibrated DDFs,
and 0.42±0.3 for HYOGA2 with globally distributed DDFs. It is clear that HYOGA2
with globally calibrated DDFs had a similar performance to the old model. Hence, theoverall time series of annual mass balance global means were similar between the oldand new models (Figure 2).

9Supplement References (for other references see main text)
Cogley JG. 2010. A more complete version of the World Glacier Inventory. Annals of
Glaciology 50: 32–38. doi: 10.3189/172756410790595859.
Dyurgerov MB, Meier MF. 1997. Mass balance of mountain and subpolar glaciers: a
new global assessment for 1961–1990. Arctic and Alpine Research 29:379–391.
Fan Y, van den Dool H. 2008. A global monthly land surface air temperature analysis
for 1948-present. Journal of Geophysical Research-Atmospheres 113: D01103.
doi: 10.1029/2007JD008470.
Haerter J, Hagemann S, Moseley C, Piani C. 2011. Climate model bias correction andthe role of timescales. Hydrology and Earth System Sciences 15: 1065–1079.
doi:10.5194/hess-15-1065-2011.
Hirabayashi Y, Kanae S, Motoya K, Masuda K, Döll P. 2008b. A 59-year (1948–2006)
global near-surface meteorological data set for land surface models. Part II:
Global snowfall estimation. Hydrological Research Letters 2: 65–69. doi:
10.3178/HRL.2.65.
Li H, Sheffield J, Wood E. 2010. Bias correction of monthly precipitation andtemperature fields from Intergovernmental Panel on Climate Change AR4models using equidistant quantile matching. Journal of Geophysical ResearchAtmospheres 115: D10101. doi:10.1029/2009JD012882.
Mitchell T, Jones P. 2005. An improved method of constructing a database of monthly
climate observations and associated high-resolution grids. International Journalof Climatology 25: 693–712. doi:10.1002/joc.1181.
Rudolf B, Rubel F. 2005. Observed Global Climate, Landolt-Börnstein (Numerical
Data and Functional Relationships), M. Hantel (eds). Springer: Berlin; 11.1111.53.
Seneviratne SI, Nicholls N, Easterling D, Goodess CM, Kanae S, Kossin J, Luo Y,
Marengo J, McInnes K, Rahimi M, Reichstein M, Sorteberg A, Vera C, Zhang X.
2012. Changes in climate extremes and their impacts on the natural physicalenvironment. In: Managing the Risks of Extreme Events and Disasters to
Advance Climate Change Adaptation. A Special Report of Working Groups I
and II of the Intergovernmental Panel on Climate Change (IPCC), Field CB,
Barros V, Stocker TF, Qin D, Dokken DJ, Ebi KL, Mastrandrea MD, Mach KJ,
Plattner G-K, Allen SK, Tignor M, Midgley PM. (eds). Cambridge University
Press, Cambridge, UK, and New York, NY, USA; 109–230.

10

Serreze MC, Walsh JE, Chapin III FS, Osterkamp T, Dyurgerov M, Romanovsky V,
Oechel WC, Morison J, Zhang T, Barry RG. 2000. Observational evidence ofrecent change in the northern high-latitude environment. Climatic Change
46:159–207.
Yatagai A, Arakawa O, Kamiguchi K, Kawamoto H, Nodzu MI, Hamada A. 2009. A
44-year daily gridded precipitation dataset for Asia based on a dense network ofrain gauges. SOLA 5: 137–140. doi:10.2151/sola.2009-035.

11

Supplement Table
Supplement Table SI. Summary of the GCMs selected in this study. The institution
and model names were taken from http://cmip-pcmdi.llnl.gov/cmip5/availability.html.
Size information was extracted from data headers.
Model

Institution

CCCma-CanESM2

Size (no. of grids)
North-South

East-West

Canadian Centre for Climate Modelling andAnalysis, Canada

64

128

CNRM-CM5

Centre National de Recherches Meteorologiques/Centre Europeen de Recherche et Formation
Avancees en Calcul Scientifique, France

128

256

GFDL-ESM2G

Geophysical Fluid Dynamics Laboratory, USA

90

144

INM-CM4

Institute for Numerical Mathematics, Russia

120

180

MIROC5

Atmosphere and Ocean Research Institute, National
Institute for Environmental Studies, and Japan
Agency for Marine-Earth Science and Technology,
Japan

128

256

MPI-ESM-LR

Max Planck Institute for Meteorology (MPI-M),
Germany

96

192

MRI-CGCM3

Meteorological Research Institute, Japan

160

320

NCAR-CCSM4

National Center for Atmospheric Research, USA

192

288

NCC-NorESM1-M

Norwegian Climate Centre, Norway

96

144

12

a) HYOGA

25
20
15
10
5

n= 46
avr.: 0.47
std.: 0.24

0
-1.0

Number of stations

Number of stations

30

30

-0.5

0.0

0.5

correlation coefficientb) HYOGA2

25
20
15
10
5

n= 71
avr.: 0.42
std.: 0.30

0
-1.0

-0.5

0.0

0.5

correlation coefficient30c) HYOGA2 with calibration

25
20
15
10
5

n= 114
avr.: 0.48
std.: 0.25

0
-1.0

1.0

Number of stations

Number of stations

Supplement Figures

30
25

0.0

0.5

1.0

d) HYOGA2 with precipitation
lapse rate

20
15
10
5

n= 114
avr.: 0.46
std.: 0.25

0
-1.0

1.0

-0.5

correlation coefficient

-0.5

0.0

0.5

correlation coefficient

1.0

Supplement Figure S1. Frequency of correlation coefficients between observed andmodeled annual mass balances of the old model HYOGA (a), HYOGA2 with globally
distributed parameters (this study; b), HYOGA2 with calibration (c), and thesensitivity test of HYOGA2 with a 10% (100m)-1 increase in precipitation withaltitude (d). Number of glaciers (n), averages (avr.,) and standard deviations (std.) ofcorrelation coefficients among all stations are given in the figures.

13