Projection of glacier mass changes under a high-emission climate scenario using the global glacier model HYOGA2

: We report a time series (1948–2100) of global-scale meltwater from mountain glaciers and ice caps (MGI) estimated by the global glacier model HYOGA2. HYOGA2 calculates the temporal fluctuation of the mass balance for 24,234 individual glaciers worldwide. It covers 90% of the total glacier area, except for glaciers in Greenland and Antarctica. HYOGA2 also accounts for regionally distributed changes in glacier area and altitude associated with glacier retreat and advance. By computation of individual glacier changes, future dissipation and glacier mass and area changes can be simulated in the model. The cumulative volume loss of water between 1948 and 2005 was estimated to be 25.9 ± 1.4 mm sea level equivalent (SLE). A future projection under a high-emission scenario demonstrated significant losses of water from MGI equivalent to 60.3 ± 7.9 mm SLE between 1948 and 2060 and 99.0 ± 14.9 mm SLE between 1948 and 2099.


INTRODUCTION
Meltwater from mountain glaciers and ice caps (MGI) is an important water resource in arid and semi-arid river basins mainly fed by mountainous water storage such as snowpack and ice (IPCC, 2007;Immerzeel et al., 2010;Kaser et al., 2010). Although measured and model-based estimations of the meltwater contribution to water resources vary widely among regions, MGI makes a significant contribution to streamflow on continental to global scales. This highlights the potential vulnerability of water resources under ongoing climate warming and subsequent shrinkage of MGI in some regions (Immerzeel et al., 2010;Kaser et al., 2010). Observations show that the present rise in runoff is a response to glacier melting in some glacier-fed basins in Alaska, the tropical Andes, and the Alps (IPCC, 2007). As the climate becomes warmer, river discharge will initially increase due to enhanced melting, until decreases in glacier volume results in reduced water yields (Braun et al., 2000;Jansson et al., 2003). A numerical future projection of terrestrial glacier mass change under a warming climate at selected large Asian rivers suggested that the disappearance of glaciers would lead to a decrease in river discharge. This could cause a loss in food production affecting as many as 60 million people (Immerzeel et al., 2010).
Early continental-to global-scale models of glacier mass change treated the glacier area distribution as steady (Hock et al., 2009;Kaser et al., 2010). However, modeling the evolution of glacier area and its altitudinal allocation, based on changes in glacier volume due to changes in atmospheric forcing is indispensable for projecting future glacier mass changes. Because the future change in climate variables simulated by General Circulation Models (GCMs) varies in different regions, especially for precipitation, estimation of the mass change in individual glaciers with a distributed glacier model is required to estimate the climate change impact on water resources in glacier-fed river basins. In particular, analysis of regime changes (i.e. increased runoff when glaciers are melting turning to decreased runoff as glaciers begin to disappear) of basin-scale glacier runoff using a model that can calculate the dynamic evolution of glacier area is crucial for future water resources assessment.
Using a regionally distributed glacier model, Radić and Hock (2011) estimated future sea level rises from the regionally distributed mass change in MGI under climate projections by several GCMs. They modeled the individual MGI using seven model parameters obtained from parameter transfer functions calibrated at 36 glaciers located in different climate and topographic regions. Possibly because of the lack of observed MGI mass changes, however, the modelestimated volume changes were not validated and only future projections were presented. Marzeion et al. (2012) calculated the global surface mass balance of glaciers in the past and future based on a complete glacier inventory, the Randorf Glacier Inventory (RGI), and monthly precipitation and temperature. For integration into internationally significant reports such as the assessment report of the Intergovernmental Panel on Climate Change (IPCC), it is necessary to assign uncertainty ranges to model-simulated future projections of glacier volume by carrying out additional numerical simulations using independent models.
A simple temperature-index-based global glacier model, HYOGA, was developed to simulate the current and future mass changes in MGI and was validated globally (Hirabayashi et al., 2010). Compared to available observations, the HYOGA simulation, calculating the mass changes of MGIs aggregated into 0.5° model grid cell at 50 m vertical resolution, reproduced reasonably well the decadal glacier mass changes and the recent (after the 1990s) acceleration of glacier mass loss. Hirabayashi et al. (2012) updated the model to estimate the mass balances of individual MGI by including a glacier inventory. They also validated the model performance at glaciers across Europe. The present study describes an updated version of the model, HYOGA2, which utilizes an extended global-scale glacier inventory, the RGI. Future changes in glacier volume are then projected using HYOGA2 forced by the climate projections of nine independent GCMs.

Global glacier model HYOGA2
HYOGA2 is an updated version of the global glacier model HYOGA (Hirabayashi et al., 2010). The first version of HYOGA2 was validated over glaciers in Europe . The current study extends the results of Hirabayashi et al. (2012) to a global scale, including the complete glacier inventory of the RGI and using globally distributed parameters obtained through a global calibration by a method similar to that suggested by Hirabayashi et al. (2010). In this section, the basic structure and updated components of the HYOGA2 model are briefly described.
The basic model structure of HYOGA2 is the same as that of HYOGA. The melt rates on snow and ice surfaces at each elevation band are calculated using a simple degreeday approach and then aggregated to estimate the mass balance for each MGI. The precipitation and surface temperature at a daily temporal resolution are necessary inputs. The subgrid-scale variation in glacier and snow mass changes is considered by dividing a model grid cell into 50 m elevation bands. The surface air temperature is assumed to decrease at a constant lapse rate of −0.65°C (100 m) −1 . Precipitation is assumed to be constant for all elevation bands in each grid cell because an approach suitable for determining precipitation lapse rates on a global scale is lacking. Although precipitation lapse rates may affect the model-estimated mass balance in a warming climate, we did not use it to avoid the incorporation of an additional uncertain parameter. Nonetheless, we note that the incorporation of adequate precipitation lapse rates is an important subject of future research. A detail discussion is given in Supplement Information S1. Further studies are expected to examine the effect of precipitation lapse rates in a future warming climate and to set it adequately on a global scale.
The glacier volume was obtained using the equation where A (km 2 ) and V (km 3 ) are, respectively, the total glacier area and volume in each glacier, and c a and γ are model parameters. For all glaciers, the model assumed fixed parameter values: γ = 1.375 from Bahr et al. (1997) and c a = 0.2055 (m 3−2γ ) from Chen and Ohmura (1990). The original model (HYOGA) assumed a single ice mass in each model grid cell, with the area determined from the corresponding total glacier area in the inventory. However, the new HYOGA2 independently calculates the mass balances of each MGI using information from the existing glacier inventory. In addition, differently from HYOGA, HYOGA2 updates the maximum glacier length L after each time step based on the new total glacier volume. This is defined following the method of Radić and Hock (2010) as where c i (= 1.7026 m 3−2γ ) and q (derived from A and L from glacier inventory) are fixed model parameters.
Because of the nonlinear relationship between A and V in Equation (1), HYOGA, which considers a single ice mass with large area, potentially overestimated the total volume of the glaciers. HYOGA2 overcomes this deficiency by separately simulating individual glaciers with relatively smaller areas. As a result, the potential overestimation of the volume of each small glacier is relatively low. More importantly, with the inclusion of independent glacier volume estimations from the completed glacier inventory, the new HYOGA2 can estimate glacier dissipation under warming climate scenarios globally.
In HYOGA2, mass changes are calculated for all the glaciers, excluding those in Greenland and Antarctica larger than 2 km 2 , to reduce the calculation cost. MGI with areas < 2 km 2 (LT2) are then calculated by an upscaling method. The LT2s, corresponding to about 10% of the total glaciers by area, are important because they are mostly located at lower altitudes, where the glacier mass change is the most sensitive to climate change (Oerlemans and Fortuin, 1992). First, the LT2s within a 0.5° grid cell were ordered according to area. Then, the median LT2 was selected as a representative LT2 for the grid. The LT2s smaller than 0.01 km 2 , which were normally distributed within a few modeled elevation bands of 50 m resolution, were not taken into consideration for the selection, since their inclusion resulted in model instability because of the small number of elevation bands. The mass change per unit area of a representative LT2 was calculated and then multiplied by the total LT2 area to obtain the total mass change of all LT2s for each 0.5° grid cell. This process assumed that all small LT2s located in the same 0.5° grid showed the same mass changes. Evaluation of the impact of standard LT2 selection on large-scale mass changes over glaciers in Europe showed that the upscaling method used in this study resulted in very similar (< 0.01%) total ice mass changes compared to methods where all small individual glaciers are calculated .
There were two calibration parameters in HYOGA2: degree-day factors (DDFs) for ice and for snow. Other calibration parameters required in HYOGA were obtained from the RGI. To obtain globally calibrated DDFs, we first aggregated the glacier information from the RGI into global 0.5° grid cells. Then, the globally distributed DDFs were obtained using a method suggested by Hirabayashi et al. (2010) at each 0.5° grid cell. Detailed descriptions of the model parameter calibration and initialization are presented in the Supplement Information S2.

Glacier inventory
The glacier inventory used in this study was version 2.0 of the Randolph Glacier Inventory (RGI, http:// www.glims.org/RGI/randolph.html) processed by Marzeion et al. (2012). The location (longitude and latitude), maximum and minimum altitude, and glacier area could be obtained from the inventory. We calculated glaciers with complete data for location and area. In total, 257,325 glaciers with a total area of 515,706 km 2 were selected from the RGI ( Figure  1). This represents an improvement compared to the original HYOGA model, which used a 0.5° gridded global glacier area dataset created by aggregating the glacier data of the older version of the World Glacier Inventory (WGI) and the Global Hydrological Data (GGHYDRO; Cogley, 2003). For details of the glacier inventories and adjustments necessary for their application in the HYOGA2 model simulations, see the Supplement Information S3.

Climate input
The retrospective simulation for a 64-year period (1948 to 2011 inclusive) of HYOGA2 was 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., 2005Hirabayashi et al., , 2008a) (details of the extension are given in Supplement Information S4). The future simulation (from 2006 to 2100) was forced by nine GCMs (Supplement Table SI) participating in the fifth phase of the Coupled Model Intercomparison Project (CMIP5), forced by the highest emission scenario (representative concentration pathway (RCP) 8.5). The extension of H08 up to 2011 and the bias correction of GCM data for 2006-2100, following the method of Watanabe et al. (2012), are briefly introduced in the Supplement Information S4 and S5.

RESULTS
The area-weighted observation-based annual mass balances and old and new model results at regional to global scale are presented in Figure 2. For the global mean, both simulations show a global acceleration of negative mass balance after the 1990s, as also reported by Dyurgerov and Meier (2005) (DM05). It should be noted, however, that the regional estimation of DM05 with a limited number of mass balance observations does not necessarily reflect the real  mass balance changes over the region. Even with similar model parameters, a significant difference between the old and new model may be due to the inclusion of additional glaciers that were missing in the old model. In North America (NA), many glaciers in arctic Canada were newly included in HYOGA2. Because of the small number of observations there, however, the calibration of DDFs was performed by assuming that the glaciers were in balance between 1948 and 1980. This caused high (closer to 0) mass balances in the model simulations because of the inclusion of more glacier information in these arctic regions (e.g., the arctic Canada in NA). In South America (SA), the overall time series did not match with the observation-based estimation. This was mainly due to the limited number of mass balance observations, in addition to the relatively lower quality of the climate input data.
The mean and standard deviation of correlation coefficients between the modeled and observed mass balance among available observation stations showed that HYOGA2 with globally calibrated DDFs had a similar performance to HYOGA (Supplement Figure S1). The root mean-square error between the results of the retrospective simulation with HYOGA2 and the available observed mass balance was 915 mm. This was close to the range of the error estimation for a global glacier model reported by Marzeion et al. (2012). Comparison of modeled and observation-based areaweighted regional averages of mass changes by Dyurgerov and Meier (2005) revealed that decadal trends of modeled mass changes on regional scales were similar between the old and new models, but the long-term mean of annual mass balances was different. This may have been due to the inclusion of new inventory. Nonetheless, the global averages of the two models were similar (Figure 2).
Even though the results of the retrospective simulations were similar to those of old models, HYOGA2 is advantageous for estimating future volume change because it distinguishes individual glaciers. In particular, HYOGA2 can estimate the retreat and dissipation of glaciers under a warming climate more adequately than the original HYOGA model. This is important for analyses of regime changes (i.e., rise in runoff turns to fall in runoff as glaciers diminish) of basin-scale glacier runoff for projecting future water resources.
A future projection of cumulative glacier mass loss and its contribution to future sea level rise based on HYOGA2 is summarized in Table I and Figure 3. A retrospective HYOGA2 simulation calculated that the global total volume loss from glaciers for 1948-2005 was 9.4 ± 0.5 × 10 3 km 3 , corresponding to a 25.9 ± 1.4 mm sea level equivalent (SLE). According to the HYOGA2 simulation forced by the outputs of nine GCMs under the most severe warming scenario (RCP 8.5), the cumulative volume loss from global glaciers becomes 21.8 ± 2.9 × 10 3 km 3 by 2060 and 35.8 ± 5.4 × 10 3 km 3 by 2099. Consequently, the total sea level rise since 1948 due to meltwater from MGI was estimated to be 60.3 ± 7.9 mm by 2060 and 99.0 ± 14.9 mm by the end of the 21 st century. These estimations are relatively smaller than that of Marzeion et al. (2012), who showed 217 ± 47 mm up to 2100 with 15 GCMs forced by the RCP8.5 scenario.
The time series of modeled mass balance showed an obvious acceleration of the negative trend since around 2000. This trend continued until the end of the simulation (Figure  3a), as a result of the assumed warming climate in the future projection. The cumulative volume loss from glaciers therefore increased rapidly, especially after the middle of the 21 st century (Figure 3b).

DISCUSSION
Degree-day factors for ice and snow may show  Table I. Future projection of the average volume loss of MGI since 1948 (10 3 km 3 ) and its sea level equivalent (mm). The value for 1948-2005 is the result of retrospective simulations with root mean-square errors of available observations multiplied by the total glacier area, excluding Greenland and Antarctica (5.17 × 10 5 km 2 ). The values in future periods are means and standard deviations among the GCMs. The ocean surface area is assumed to be 362 × 10 6 km 2 .

Period since 1948
Volume loss since 1948 (10 3 km 3 ) Sea level equivalent (mm) -2005 9.4 ± 0.5 25.9 ± 1.4 -2060 21.8 ± 2.9 60.3 ± 7.9 -2099 35.8 ± 5.4 99.0 ± 14.9 significant spatial variability (Zhang et al., 2006), greatly affecting the accuracy of snow-and ice-melt models. As in the old HYOGA model, DDFs were calibrated until they reproduced similar values for the grid-cell-specific longterm average of total glacier mass balance observed in neighboring grid cells (Hirabayashi et al., 2010). Because the DDFs were the only calibration parameters of HYOGA2, the calibrated DDFs might have compensated for other potential sources of errors in the model simulations. For example, because of the limitations of the global climate input data, the model simulation might not be accurate at glaciers located in poorly gauged regions, such as High Mountain Asia. Glacier characteristics, which were not adequately represented in the model (e.g., slope aspects and effects of debris cover), were also potential error sources that might have partly been compensated for in the DDF calibration.
In addition, debris-covered glaciers are common in many mountain ranges around the world. In these glaciers, the debris cover in the ablation zones accelerates (if the debris is thin) or suppresses (if it is thick) the melting of sub-debris ice relative to that of other snow and ice (Nakawo andYoung 1981, Mattson et al. 1993). This then affects the DDFs (Zhang et al., 2006). However, the effect of the debris cover on DDFs could not be considered in HYOGA2 and other global glacier models (Radić and Hock, 2011;Marzeion et al., 2012) because of the poor understanding of the largescale spatial distribution of debris thickness and properties. One potential future method to obtain these debris-related properties is to use remote sensing data to relate the debris thickness and thermal resistance of the debris layer (Zhang et al., 2011).
Improper model initialization is also known to cause biases in model simulation. Because we assumed that the glacier area was provided in the RGI as an initial area in 1948, the glacier area may often be underestimated for glaciers where melting was already obvious in the 20 th century, if RGI data were recorded at the end of the 20 th century. With respect to future projection in the 21 st century (Figure 3), we assume that this effect is small compared to the computed future mass loss. The upscaling method for LT2s, which assumes the same mass changes, may also cause potential errors.
Finally, in addition to glacier model-related sources of errors, other potential sources include the model input data and processing. Future GCM simulations that are used to force glacier models have significant biases and require correction. Future studies should therefore also focus on assigning uncertainty ranges to model simulations because of uncertainties associated with the selection of bias correction methods for GCM outputs.

SUMMARY
A global glacier model, HYOGA2, which includes a completed individual glacier inventory and uses globally distributed parameters, was developed and used to estimate the global glacier mass changes in the past and future. The performance of the retrospective HYOGA2 simulation was similar to that of HYOGA, showing reasonably accurate temporal variations compared to observations at both global and regional scales. Because of the inclusion of the complete global-scale glacier inventory, the updated model could represent large-scale glacier mass losses in the future better than the old model. Dissipation of small glaciers under a warming climate is one of the key factors to consider in the projection of future changes in glacier meltwater volume. By calculating the mass changes in individual glaciers, HYOGA2 can simulate the timing of MGI dissipation and the associated shortage in water resources under future climate scenarios. Future model simulations using nine GCMs forced by the highest emission scenario (RCP 8.5) showed significant decreases in water from MGI, with losses amounting to 60.3 ± 7.9 mm SLE between 1948 and 2060 and 99.0 ± 14.9 mm SLE between 1948 and 2099.

ACKNOWLEDGMENTS
This research was supported by the Funding Program for next-generation world-leading researchers, the Japan Society for the Promotion of Science, and the CREST project of the Japan Science and Technology Agency.

SUPPLEMENTS
Detail explanations of model sensitivity to precipitation lapse rate (Supplement Information S1), model calibration and initialization (Supplement Information S2), glacier inventory (Supplement Information S3), extension of observation-based climate data (Supplement Information S4), bias correction of the GCM data (Supplement Information S5) and model validation (Supplement Information S6 and Supplement Figure S1) are given in the supplementary materials. Summary of the GCMs used in this study is given in the Supplement Table SI.