Past and future sea-level change from the surface  mass balance of glaciers

<jats:p>Abstract. We present estimates of sea-level change caused by the global surface mass balance of glaciers, based on the reconstruction and projection of the surface mass balance of all the individual glaciers of the world, excluding the ice sheets in Greenland and Antarctica. The model is validated using a leave-one-glacier-out cross-validation scheme against 3997 observed surface mass balances of 255 glaciers, and against 756 geodetically observed, temporally integrated volume and surface area changes of 341 glaciers. When forced with observed monthly precipitation and temperature data, the glaciers of the world are reconstructed to have lost mass corresponding to 114 ± 5 mm sea-level equivalent (SLE) between 1902 and 2009. Using projected temperature and precipitation anomalies from 15 coupled general circulation models from the Coupled Model Intercomparison Project phase 5 (CMIP5) ensemble, they are projected to lose an additional 148 ± 35 mm SLE (scenario RCP26), 166 ± 42 mm SLE (scenario RCP45), 175 ± 40 mm SLE (scenario RCP60), or 217 ± 47 mm SLE (scenario RCP85) during the 21st century. Based on the extended RCP scenarios, glaciers are projected to approach a new equilibrium towards the end of the 23rd century, after having lost either 248 ± 66 mm SLE (scenario RCP26), 313 ± 50 mm SLE (scenario RCP45), or 424 ± 46 mm SLE (scenario RCP85). Up until approximately 2100, ensemble uncertainty within each scenario is the biggest source of uncertainty for the future glacier mass loss; after that, the difference between the scenarios takes over as the biggest source of uncertainty. Ice mass loss rates are projected to peak 2040 ∼ 2050 (RCP26), 2050 ∼ 2060 (RCP45), 2070 ∼ 2090 (RCP60), or 2070 ∼ 2100 (RCP85).                     </jats:p>

section_content: By temporally integrating the surface mass balance over long periods of time, fluctuations in glacier geometries allow people to perceive slow changes of the climate system, which otherwise would be overwhelmed in human perception by short-term variability. Because of this property, shrinking glaciers around the world have become poster children of climate change. 

But impacts of glacier change -whether growing or shrinking -go far beyond this sentimental aspect: by changing the seasonality of runoff, glaciers are important regulators of water availability in many regions of the world (Kaser et al., 2010; Huss, 2011; Immerzeel et al., 2012). Retreating glaciers also lead to increased geohazards, e. g. from destabilized slopes and lakes dammed behind unstable, icecored moraines (see Richardson and Reynolds, 2000, for an overview). Finally, even though the ice mass stored in glaciers seems negligible compared to the Greenland and Antarctic ice shields, glaciers 1 have contributed significantly to sea-level rise in the past (Cogley, 2009; Hock et al., 2009; Hirabayashi et al., 2010; Leclercq et al., 2011), and probably have been the biggest single source of observed sea-level rise since 1900 (Lemke et al., 2007). 

Improving knowledge of how glaciers have been and will be changing when subjected to climate change, both natural and anthropogenic, is therefore a pressing task. The main obstacle to progress is a severe undersampling problem: direct glaciological measurements, e. g. of surface mass balances, have been performed on ∼ 300 glaciers world wide. The number of glaciers on which these types of measurements have been carried out for time periods longer than 30 yr, i. e. over periods that potentially allow for the detection of a climate change signal, is one order of magnitude smaller. Length variations of glaciers have been observed for substantially longer periods of time, and also for hundreds of glaciers (Oerlemans, 1994 (Oerlemans,, 2005)), but are much more difficult to understand, as large glacier length fluctuations may arise from intrinsic climate variability (Roe and O'Neal, 2009; Roe, 2011), and ice dynamics leads to a complicated relation-ship between surface mass balance and length variations (e. g. Jarosch and Marzeion, 2012). Data obtained by remote sensing (e. g. gravimetric assessments of ice mass change, or volume change estimates obtained by differencing digital elevation models) may cover a greater number of glaciers, but are available only for short and recent periods (Gardner et al., 2011; Moholdt et al., 2012; Jacob et al., 2012). But even if these difficulties were overcome, the undersampling problem would remain: the estimated total number of individual glaciers in the world is ∼ 200 000 (Radić and Hock, 2010; Arendt et al., 2012), and even the strongest effort in improving data coverage could not lead to a substantial improvement of the situation within the next decades (see Fig. 1 for an overview of the total distribution of glaciers in the world, compared to those glaciers on which measurements exist). 

Besides the practical difficulties of modeling vast numbers of glaciers, the undersampling problem imposes strong limits on (i) the number of model parameters that can be determined empirically, and (ii) the reliability and representativeness of these parameter estimates. Therefore, some recent efforts have circumvented the problem of direct modeling, e. g. by simply extrapolating observed mass loss rates and observed mass loss trends into the future (Meier et al., 2007), by prescribing future surface mass balance and discharge rates and testing the plausibility of meeting the demands that follow from the prescribed rates (Pfeffer et al., 2008), or by assuming either constant or constantly declining accumulation area ratios in the future (Bahr et al., 2009). Raper and Braithwaite (2006) model future mass balances of glaciers by first deriving statistical characteristics of the glacier distributions within grid cells of 1 × 1 degree, then using climate data to derive mass balance profiles for each grid cell, but ultimately have to rely on the extrapolation of the results from seven geographically and climatically limited regions with enough data coverage to the rest of the world.
section_name: Introduction

section_content: Of all published studies, Radić and Hock (2011) employ the most complex surface mass balance model, individually for each known glacier (which at the time of publication implied that only about half of the world's glacierized area apart from the ice sheets was modeled directly). Their model parameters are determined from a rather small number of glaciers with measured vertical mass balance profiles, and because they ultimately tune one parameter to fit their model results to the observed, regionally integrated mass balances of Dyurgerov and Meier (2005), an independent validation of the model set up during the 20th century is not performed. 

In summary, even though impacts of glacier change affect people much more directly than changes of the Greenland and Antarctic ice sheets, and even though glaciers will potentially contribute more strongly to sea level rise within the 21st century than the ice sheets (Meehl et al., 2007), very few studies have projected future change of the world's glaciers. To our knowledge, no global projection has been performed using a model that was validated independently against observed glacier changes of the 20th century. But an independent validation of model systems that are used to project the future behavior of (parts of) the climate system is paramount for the trustworthiness of the projections (see e. g. Randall et al., 2007). This is particularly true if the knowledge about the system to be modeled is limited, as is the case with the world's glaciers. 

Here, we will present reconstructions and projections of glacier change based on a model system that tries to make best use of the relatively few measurement points that are available for model validation, leading to a statistically robust assessment of the model's skill and errors. Additionally, the model system allows for seamless simulation of past and future glacier changes, in order to ensure that the future skill and error properties of the model projection can reliably be determined from the past. 

In Sect. 2, we present the mass balance model, which includes a simple representation of glacier geometry change in response to climatically forced volume changes, and describe the sources and estimations of the necessary model parameters and variables. Then, the model is first applied for the 20th century, and two independent validations are performed, where we deliberately put the emphasis on the validation of the model's results, rather than on the accuracy of the parameter estimates: in Sect. 3, a leave-one-glacier-out cross-validation of the modeled surface mass balance of 255 glaciers with measured surface mass balances (blue dots in Fig. 1 ) is presented. Based on the results from this crossvalidation, the model error is propagated through the entire model in order to obtain uncertainty estimates for each of the modeled variables (Sect. 4). Then, modeled volume and surface area changes, as well as their uncertainty estimates, are validated again using geodetically measured volume and surface area changes of 341 glaciers (green rings in Fig. 1 ) in Sect. 5. We describe the data sets used for forcing the model in Sect. 6. Finally, the results are presented in Sect. 7, and discussed in the context of other reconstructions and projections in Sect. 8.
section_name: B. Marzeion et al.: Sea-level change from glaciers 1297

section_content: For each individual glacier, we calculated the annual specific surface mass balance B as

where P solid i is the area mean monthly solid precipitation onto the glacier surface (see Sect. 2. 1. 1), µ * is the temperature sensitivity of the glacier (see Sect. 2. 1. 9), T terminus i is the monthly mean air temperature at the location and elevation of the glacier's terminus (see Sect. 2. 1. 2), T melt is the monthly mean air temperature above which ice melt is assumed to occur (see Sect. 2. 2. 5), and β * is a bias correction (see Sect. 2. 1. 10). We thus do not attempt to capture the full energy balance at the ice surface, but rely on air temperature as a proxy for the energy available for melt (Ohmura, 2001; Hock, 2003; Sicart et al., 2008).
section_name: Mass balance model

section_name: Glacier-specific model parameters and variables

section_content: The area mean monthly solid precipitation onto the glacier surface P solid i is estimated as

where a is a precipitation correction factor (see Sect. 2. 2. 2); P CRUclim i is the monthly total climatological precipitation taken from the grid point of the CRU CL 2. 0 dataset (10 spatial resolution, New et al., 2002) closest to the glacier; P anom i is the monthly total precipitation anomaly taken either from the closest grid point of the CRU TS 3. 0 dataset (0. 5 • spatial resolution, Mitchell and Jones, 2005) (for the case of applying observed climate variability and change), or from the closest grid point of the climate model (for the case of applying modeled climate variability and change, see Sect. 6. 1 for an overview over the applied models and scenarios); γ precip is a precipitation lapse rate (see Sect. 2. 2. 3); z mean is the mean elevation of the glacier (see Sect. 2. 1. 7); z CRUclim is the elevation of the grid point of the CRU CL 2. 0 data set; and f solid is the fraction of solid precipitation. 

The fraction of solid precipitation is based on the monthly mean temperature as where T prec solid is the temperature below which precipitation is assumed to be solid (see Sect. 2. 2. 4), γ temp is a temperature lapse rate (see Sect. 2. 1. 8), z max is the maximum elevation of the glacier (see Sect. 2. 1. 6), and z terminus is the terminus elevation of the glacier (see Sect. 2. 1. 5). This implies that all precipitation falling onto the glacier surface is assumed to be solid if the monthly mean temperature at the terminus elevation of the glacier is below T prec solid, that all precipitation is assumed to be liquid when the monthly mean temperature at the maximum elevation of the glacier is above T prec solid, and that the fraction of solid precipitation decreases linearly with temperature between these two points.
section_name: Precipitation

section_content: The monthly mean air temperature at the location and elevation of the terminus of the glacier is estimated as

where T CRUclim i is the monthly mean climatological temperature taken from the grid point of the CRU CL 2. 0 dataset closest to the glacier; and T anom i is the monthly mean temperature anomaly taken either from the closest grid point of the CRU TS 3. 0 dataset (for the case of applying observed climate variability and change), or from the closest grid point of the climate model (for the case of applying modeled climate variability and change).
section_name: Temperature

section_content: The measured surface area A measured of a glacier is determined by integrating the glacier outlines from the Randolph Glacier Inventory (Arendt et al., 2012) (RGI) version 1, which is a data set of glacier outlines which combines previously existing outlines from different sources (including the Global Land Ice Measurements from Space (GLIMS) initiative data) with new data from various contributors. It is the first globally complete glacier inventory. 

For those regions in which individual glaciers are not separated in the RGI, we drape the RGI glacier outlines over the version 2 of the ASTER global digital elevation model (GDEM) and use a watershed algorithm (Ehlschlaeger, 1989) to separate the outlines into individual glaciers. After identifying drainage basins on the whole GDEM which are at least 9 km 2 in size 2, we crop these basins with the RGI glacier outlines to delineate individual glaciers. This is possible because, at least as a zero-order approximation, the same physics underlie the identification of drainage basins and individual glaciers, even though the non-linearity of ice rheology is quite complex. If the topography did not allow for a clear separation into drainage basins, then the outline was examined manually, and, if appropriate, treated as an ice cap (see Sect. 2. 2. 1). Overall, 29 ice caps were identified this way. Since no exact dates of the determination of glacier outlines are given in the RGI, we estimate approximate years of area measurement from the technical document accompanying the RGI data set, and account for the additional uncertainty in the timing of the area measurement by adding uncertainty to the area measurement (see Sect. 4). 

The surface area change dA of the glacier during each mass balance year (running from October to September in the Northern Hemisphere, and April to March in the Southern Hemisphere) is modeled as

where τ A is a relaxation time scale (see Sect. 2. 1. 12), V (t +1) is the glacier's volume at the end of the mass balance year (see Sect. 2. 1. 4), c A and γ are scaling parameters (see Sect. 2. 2. 1), and A(t) is the surface area of the glacier at the end of the preceding mass balance year. This means that an equilibrium surface area (i. e., V (t+1) c A 1/γ in Eq. 5) is estimated corresponding to the volume of the glacier modeled (Bahr et al., 1997; Bahr, 1997). But the surface area of the glacier does not take this equilibrium value instantly, instead it is linearly relaxed towards it from its current surface area. 

Since the relaxation time scale τ A introduces memory of past changes into the model, it is not possible to integrate the model backwards in time to determine the evolution of the glacier before the year of surface area measurement. For this reason, the glacier's surface area A start at the beginning of the model integration (i. e. 1901 for the forcing with observed climate variability and change, and 1850 for most cases of modeled climate variability and change) is estimated by iteratively seeking that surface area in the starting year of the integration that will result in the measured surface area in the year of the measurement. The iteration is deemed successful when the modeled surface area is within 0. 1 % of the measured surface area during the year of the measurement; the iterative process is broken off after 100 iterations if unsuccessful (see Sect. 6. 2. 2 how these glaciers are treated).
section_name: Area

section_content: The glacier volume V start in the starting year of the model integration is estimated following volume-area scaling (Bahr et al., 1997; Bahr, 1997) as

The volume change dV during each mass balance year is determined as

where A(t) and B(t) are the surface area of the glacier at the start of the mass balance year and the specific mass balance during the mass balance year (see Eq. 1), respectively, and the ice density is assumed to be ρ = 900 kg m -3. We assume a linear increase of the terminus elevation z terminus with decreasing glacier length L,

where L is the glacier's length, z measured terminus is the minimum elevation of the glacier in the year of the surface area measurement, and L 0 is the length of the glacier in the year of the surface area measurement. 

z measured terminus is taken as the ASTER GDEM minimum elevation within an individual RGI glacier outline. 

At the start of the integration,

where c L and q are scaling parameters (see Sect. 2. 2. 1). During the model integration, length changes dL during each mass balance year are estimated as

where τ L is a relaxation time scale (see Sect. 2. 1. 11), and L(t) is the glacier's length at the start of the mass balance year.
section_name: Volume

section_content: The glacier's maximum elevation z max is held constant. It is determined by finding the ASTER GDEM maximum elevation within an individual RGI glacier outline.
section_name: Maximum elevation

section_content: Here, we approximate the area-mean elevation z mean as the mean of z max and z terminus. This approximation is made in order to reflect the effects of the temporal variability of a glaciers geometry on its mean elevation in a simple way.
section_name: Mean elevation

section_content: The temperature lapse rate γ temp is estimated at each glacier location by regressing temperature of 3 × 3 CRU CL 2. 0 grid points around the location of the glacier onto z CRUclim. The correlation between temperature and elevation is very high (typically > 0. 95) and above the 95 % confidence interval for all glaciers.
section_name: Temperature lapse rate

section_content: In a first step, we estimate the temperature sensitivity µ * for all glaciers with available mass balance measurements, the data being obtained from an updated version of Cogley (2009). There is a global total of 255 glaciers with mass balance records that have all the metadata needed for the parameter estimation, are covered by the CRU TS 3. 0 and CRU CL 2. 0 data sets, are indicated to be reliable by the status flag of the data set, and have at least two annual mass balance measurements. The locations of these glaciers are shown in Fig. 1 as blue dots. 

For each of these glaciers, temperature sensitivities µ(t) are estimated by requiring that We then apply Eq. ( 1 ) for each of the 255 glaciers for each of the µ(t), and for each glacier determine the year t * in which Eq. ( 12 ) is minimal

Here, B(t) modeled is the mean of the modeled mass balances during the years of mass balance measurements, B measured the mean of the observed mass balances, and β(t) is therefore the bias of the modeled mass balances. For each of the 255 glaciers with mass balance measurements, we thus determine µ * = µ(t * ) that produces the smallest possible bias β * = β(t * ). Therefore, µ(t) is not to be understood as a timevarying temperature sensitivity. t here can rather be understood as an index to different climate conditions that are realistic at the site of the glacier, and t * as the value of the index that produces the most favorable model results. The reasoning behind this procedure is that it should be possible to identify climatologies of precipitation and temperature with which the glacier is in equilibrium. Note that this does not imply that a glacier was in equilibrium with climate around the year t *, but rather that it would have been in equilibrium around the year t * if it then had had the geometry it had during the measurement of mass balances. According to our model, the vast majority of glaciers had larger extents (and thus lower termini) around t * than during the measurement period, and thus negative mean mass balances. 

For all other glaciers (red dots in Fig. 1 ), we interpolate t * from the ten closest glaciers with mass balance measurements, weighting inversely with distance, and then determine µ * by requiring that Figure 3 shows the benefit of spatially interpolating t * rather than µ * : the values of µ * obtained by interpolating t * as described above (panel b) are much more accurate than if µ * is directly interpolated within each region (panel a).
section_name: Temperature sensitivity

section_content: Since the minimal bias β(t * ) is non-negligible for some glaciers with mass balance measurements (implying that within the period of CRU TS 3. 0 data availability, we cannot identify climatologies of precipitation and temperature with which the glacier is in equilibrium), we introduce a bias correction for all other glaciers by interpolating β * from the ten closest glaciers with mass balance measurements, weighting inversely with distance.
section_name: Bias correction

section_content: The time scale of a glacier's length response to volume changes τ L is estimated as

following roughly the scaling of Jóhannesson et al. (1989). This implies that smaller glaciers react faster, and glaciers with higher mass turnover react faster.
section_name: Time scale of glacier length

section_content: Since the flow of ice of a glacier is mainly along a valley, we assume that changes in a glaciers volume are translated instantaneously into changes of its width W, following the slope of the valley sides, and estimate

2. 2 Global model parameters
section_name: Time scale of surface area

section_content: The scaling parameters for relating the equilibrium values of volume, area, and length of a glacier at equilibrium are taken from the literature as γ = 1. 375 (Bahr et al., 1997), c A = 0. 0340 km 3-2γ (Bahr, 1997), q = 2. 2 (Bahr et al., 1997), and c L = 0. 0180 km 3-q (Radić et al., 2008) for glaciers. 

For ice caps, γ = 1. 25, c A = 0. 0538 km 3-2γ, q = 2. 5, and c L = 0. 2252 km 3-q follow from assuming a radially symmetric parabolic cross section following Cuffey and Paterson (2010). 

The remaining four global parameters are optimized within their physically meaningful ranges based on the following objectives: the cross-validation results of the model (see Sect. 3) should indicate (i) that the mass balance model has a negligible global mean bias, (ii) that the variance of the measured mass balances is well captured by the modeled mass balances, and (iii) that there is no temporal trend in the error of the mass balance model. Additionally, the correlation between modeled and measured mass balances, and the model's skill score should be as high as possible, and the root mean square error as low as possible.
section_name: Volume-area and volume-length scaling parameters

section_content: There is evidence that precipitation in the CRU data sets underestimates precipitation onto glaciers (from observations, Giesen and Oerlemans, 2012 estimate a global median factor (corresponding to a in Eq. 2) of 2. 55; by optimizing a mass balance model similar to the one presented here, Marzeion and Nesje, 2012 find a mean factor of 2. 1 in Central Europe and Scandinavia). We therefore repeat the entire crossvalidation procedure (described in Sect. 3) using 1 ≤ a ≤ 3. Figure 4 shows the results of this parameter sensitivity study: while the mean correlation between modeled and observed mass balances has a maximum for a ≈ 1. 3 (panel a), and a maximum model skill is found for a ≈ 1. 9 (panel b); the modeled mass balances have a too low variability compared to observations (panel e) and a spurious trend (panel f) for a < 2. 5. For this reason, and because mean model bias and mean root mean square error (rmse) do not depend strongly on a, we set a = 2. 5.
section_name: Precipitation correction

section_content: While mean model skill shows a weak decrease with increases of the precipitation lapse rate γ precip (Fig. 5, panel b), the mean rmse tends to decrease with increasing γ precip (panel c). Most importantly, the model has a non-zero mean bias for γ precip 2 %/100 m and γ precip 4 %/100 m. We therefore set γ precip = 3 %/100 m.
section_name: Precipitation lapse rate

section_content: The temperature below which precipitation is assumed to be solid should be close to 0 • C, but probably positive. We vary 0 • C ≤ T prec solid ≤ 5 • C, and show the results in Fig. 6. The only strong dependence of model performance on T prec solid is in the mean correlation between modeled and observed mass balances, which increases for T prec solid > 2 • C (panel a). We set T prec solid = 3 • C since the mean model bias is smallest here (panel d).
section_name: Temperature threshold for solid precipitation

section_content: The monthly mean surface air temperature above which melt is assumed to occur at the glacier terminus T melt can be assumed to be close to zero. We vary -2 • C ≤ T melt ≤ 1. 5 • C (it is problematic to set T melt > 1. 5 • C, since there are glaciers that do not experience T (t * ) terminus i,clim > 1. 5 • C, implying that the estimation of their temperature sensitivity becomes impossible, see Eq. 13). Figure 7 shows that the mean correlation between modeled and observed mass balances tends to increase with increasing T melt (panel a). Also, the variance of the modeled mass balances becomes more realistic (panel e) and a spurious trend in the modeled mass balances decreases with increasing T melt (panel f). Since the model skill drops strongly for T melt > 1 • C (panel b), we set T melt = 1 • C. One could expect that a negative T melt leads to best performance, because also in months with a monthly mean temperature below freezing, melting can occur because of diurnal and intra-monthly variability. Note that the positive value applied for T melt here does not imply that ice fails to melt any given month with temperatures above freezing and below T melt, but only that the annually summed monthly mass balance, which is evaluated here, is captured best if T melt is positive.
section_name: Temperature threshold for melt

section_content: We perform a leave-one-glacier-out cross-validation (Michaelsen, 1987; Hofer et al., 2010) of the entire modeling procedure, i. e. for each of the 255 glaciers with measured mass balances, we reconstruct its mass balance for the years of measured mass balance under the assumption of not having any information besides location, surface area (and year of surface area measurement), and elevation range of that glacier, and gain a total set of 3997 pairs of annual modeled and measured mass balances, each of the modeled mass balances being derived independently of its measured counterpart. Table 1 gives a summary of the model's performance measured by the mean root mean square error (rmse), model bias, correlation between observed and modeled mass balances, and skill score (Wilks, 2006; Marzeion et al., 2012) 3, Fig. 8 shows modeled versus observed mass balances, and Fig. 9 shows the distributions of model error and model bias. 

The performance on the global scale give confidence that the model can be expected to reconstruct annual mass balances of unmeasured glaciers to a reasonable accuracy. On the regional level, however, large discrepancies in performance become apparent. Most importantly, the mean bias of some regions does not seem to be small (see Fig. 10 ), which indicates that accumulating modeled mass balances forward in time may be problematic. But this issue is most likely the result of too small sample sizes. For only one of the regions (Scandinavia) does a Kolmogorov-Smirnov goodness-of-fit test reject (at the 95 % confidence level) the null hypothesis that the bias values within the region are drawn from a normal distribution with zero mean and the standard deviation equal to the standard deviation of the model biases in that region. In the case of Scandinavia, it is not the mean of the biases (40 mm w. e. yr -1 ), but kurtosis and a slight skewness of the distribution that are detected by the test. There is no spurious trend in the modeled mass balances (Fig. 11 ). 

The model error does not depend on glacier surface area or the number of mass balance measurements available for the cross-validated glacier (Fig. 12, panels a and b). This indicates that the model is robust, and any potential unrepresentativeness of the sampled glaciers does not affect the model's performance. Perhaps not surprisingly, there is a dependence of the root mean square error and magnitude of the bias on the remoteness of the glacier, i. e. on the mean distance to the ten closest measured glaciers that were used to determine t * (note that this indicates that, particularly for remote glaciers, the cross-validation probably underestimates the model's performance, as it artificially increases the remoteness of the glaciers in the region that is being tested). However, that dependence (particularly of the bias) is surprisingly weak (Fig. 12, panels c and d ). 

All these metrics of the model's performance proved to be relatively insensitive to parameter choices such as changing the number of closest glaciers used to determine t *, weighting the interpolation of t * by the inverse of β * instead of the distance, or weighting by a combination of distance and β *, etc. But increasing the number of required mass balance measurements for a glacier to be included in the network of glaciers used to interpolate t * (see Sect. 2. 1. 9), and thereby reducing the number of glaciers in that network, substantially increased the root mean square error of the modeled mass balances.
section_name: Cross-validation of mass balance model

section_content: Uncertainty enters the model in several places: (i) uncertainty of the mass balance model itself, (ii) uncertainty of the forcing of the mass balance model, (iii) uncertainty in the surface area measurement and in the measured maximum and terminus elevations, (iv) uncertainty in the scaling relationships relating glacier surface area to glacier volume and length, and (v) uncertainty in the representation of dynamic glacier response to volume changes, i. e. in the response time scales.
section_name: Treatment of uncertainty

section_content: In principle, the errors associated with uncertainty sources (i) and (ii) are determined directly by the cross-validation routine. However, the periods during which mass balance measurements are available are short (mean length of observation ≈ 15 yr) compared to the period over which the model is applied (≈ 100 yr for reconstructions with CRU data, and > 200 yr for projections with the GCMs). Changes in a glacier's terminus elevation, which affect the temperature that it experiences, act on longer time scales, so model uncertainties related to uncertainty in terminus elevation is probably underestimated by the cross-validation. We therefore treat this source of uncertainty separately, by first translating the uncertainty in glacier length L into uncertainty in terminus elevation z terminus, and subsequently into uncertainty in temperature T terminus, by using the temperature lapse rate γ temp determined for that glacier as described above. For each of the n months with temperatures above T melt at z terminus, this translates into an uncertainty of µ *

• (T terminus ), where (T terminus ) is the uncertainty in T terminus. For any single year, the total uncertainty of the modeled mass balance (B modeled ) then is given by

Note that while this implies that the uncertainty of the modeled mass balance generally grows over time both before and after the date of measurement of the glacier's surface area, the fact that climate variability changes the number of months n with temperatures above T melt at z terminus will lead to interannual variability in (B modeled ), and a trend in n (e. g. following anthropogenic warming) will lead to a trend in (B modeled ).
section_name: Mass balance

section_content: In the year of the surface area measurement, uncertainty of the volume, (V ), is estimated as

where VA is the relative uncertainty assigned to the volumearea scaling, estimated to 40 % 4. Before and after the year of the surface area measurement, all errors are propagated through the entire model, assuming a normal distribution, and no temporal correlation of the model errors (as indicated by the results of the cross-validation, see Fig. 9 ).
section_name: Volume

section_content: The uncertainty of the surface area measurement is generally assumed to be small. However, because no exact date of surface area measurement is given with the RGI data, we estimated dates from the description of the RGI data set, and set the uncertainty of the measured surface area to 5 %, in order to include a potential error caused by mis-dating the surface area measurement. Before and after the year of the surface area measurement, all errors are propagated through the entire model, under the same assumptions as above (Sect. 4. 2).
section_name: Surface area

section_content: In the year of the surface area measurement, the uncertainty of the glacier length (L) is estimated as

where VL is the relative uncertainty assigned to the volumelength scaling, estimated to be 100 %. Before and after the year of the surface area measurement, all errors are propagated through the entire model, under the same assumptions as above (Sect. 4. 2).
section_name: Length

section_content: The relative uncertainty of the response time scale of a glacier's surface area and length to changes in volume is estimated high at 500 %, following the analysis of explicitly modeled response times of the glaciers, integrating an ice dynamics model of a glacierized mountain range over > 1000 yr (Jarosch and Marzeion, 2012). Even so, the uncertainty that enters our model through the time scales of glacier response is small compared to those entering through the mass balance, and the volume-area and volume-length scaling.
section_name: Time scales of glacier response

section_content: Uncertainty in the measured maximum and terminus elevations was ignored, as it is negligible compared to the other sources of uncertainty. More importantly, there is probably unquantifiable uncertainty due to sampling issues: generally speaking, the regional density of mass-balance observations is well correlated with the density of weather observations (> 50 % of the glaciers with mass balance measurements are situated in the regions Scandinavia, Western Canada and US, and Central Europe, which all have a relatively high density of meteorological stations). Therefore, we can assume that the uncertainty in the forcing of the model is on average lower at the glaciers with mass balance measurements than at glacier sites in general. Similarly, most of the mass balance measurements were made in the second half of the 20th century, when the density of weather stations was higher than in the first half of the 20th century. Therefore, it is likely that the cross-validation underestimates the model error, but it is not practicable to quantify these sources of error. Note however that the validation of the integrated modeled volume and surface area changes (see Sect. 5) indicates that these sources of error are probably small.
section_name: Unquantified errors

section_content: In order to validate the modeled, temporally integrated changes of glacier volume and surface area, as well as the propagated model errors, we model each of the glaciers from Cogley (2009) for which geodetic volume change measurements exist, for which all necessary metadata are available, and which are covered by CRU data (see Sect. 6. 2 for how data gaps are treated in general). Altogether, there are 341 such glaciers (green markers in Fig. 1 ) with 756 geodetically measured volume changes. We model each of these volume change measurements twice: once using the surface area measurement of the glacier from the starting date of the volume change measurement, i. e. running the model in forecast mode; and once using the surface area measurement of the glacier from the ending date of the volume change measurement, i. e. running the model in hindcast mode. Figure 13 shows the results for one geodetic volume change measurement from Hintereisferner as an example. For each of the glaciers, and each of the geodetic volume change measurements, we then determine the difference between the modeled volume change and observed volume change, and between modeled surface area change and observed surface area change. We divide these differences by the propagated model error in order to obtain the relative model error. Over all the modeled volume changes and surface area changes, the relative error should have a mean of zero and a standard deviation of one if the modeled volume and surface area changes, and the propagated model errors, are correct. 

Figure 14 shows the distributions of the relative errors. For non-calving glaciers, the mean (0. 14) and median (0. 12) of the relative volume error are reasonably small. Positive values indicate that the glaciers were losing more volume (or surface area) than predicted by the model, or gaining less volume (or surface area) than predicted by the model. Since the mean bias of the modeled surface mass balance is very close to zero (see Sect. 3), this indicates that internal or subglacial processes are responsible for this mass loss. The mean of the relative volume error can be brought very close to zero by artificially subtracting 70 mm w. e. from the annual mass balance of each glacier every year, giving an indication of the magnitude of the volume loss through internal or basal melt processes. 

The 15th and 85th percentiles of the relative volume error lie within -1 to 1, indicating that the propagated model error is slightly larger than justified by the comparison with the geodetic volume changes. While the 2nd percentile is greater than -2, a skewness of the distribution (caused exclusively by glaciers in the Southern Andes and Alaska) causes the 98th percentile to be greater than 2. For calving glaciers, the skewness is even more pronounced, the mean (0. 74) and median (0. 35) values of the relative volume error are more positive, and the distribution of the relative volume error is wider. This was to be expected, since the model ignores solid ice discharge, which contributes to the volume changes of these glaciers. 

The model slightly overestimates surface area losses. The mean (-0. 98 for non-calving, and -0. 07 for calving glaciers) and median (-0. 11 for non-calving, and -0. 10 for calving glaciers) relative errors are negative. The distribution for non-calving glaciers has long tails, which causes the 2nd percentile to be smaller than -2, and the 98th percentile to be larger than 2. 

All in all, this implies that the modeled surface area changes are less reliable than modeled volume changes. It also implies that the propagated errors of the model for both volume and surface area changes overestimate the model uncertainty at one standard error, and underestimate the model uncertainty at two standard errors -i. e. there is excess kurtosis in the error distribution. 

Because of the slight underestimation of volume losses, there is a weak but significant correlation between the relative volume change error and the length of the time span covered by the geodetic volume change measurement for both calving and non-calving glaciers (Fig. 15 ), but this correlation disappears in hindcast mode for non-calving glaciers. Similarly, the overestimation of surface area loss leads to a weak (but in the case of non-calving glaciers significant) anti-correlation of the relative surface area change error with the time span covered by the geodetic volume change measurement. However, as already indicated by the results of the cross validation (Sect. 3), Fig. 16 shows that there is no correlation between volume change or surface area change errors and the glacier surface area (implying that any potential unrepresentativeness of the sampled glaciers does not matter). Also the remoteness of the glacier has only weak influence on the model error -the only significant correlation is between relative surface area change error and mean distance to the 10 closest sampled glaciers for calving glaciers. 

6 Forcing data and treatment of data gaps
section_name: Validation of 20th century model results

section_content: We apply the model for all individual glaciers, first using CRU precipitation and temperature, and then using output from 15 CMIP5 models as driving data set. Table 4 gives an overview of the applied climate models and scenarios (Representative Concentration Pathways, RCPs), and the respective periods. The RCPs are named after the likely anthropogenic radiative forcing of the atmosphere they represent in the year 2100, ranging from 2. 6 to 8. 5 W m 2 (see van Vuuren et al., 2011, for an overview). 

For each of the models, and each of the scenarios, anomalies of precipitation and temperature were calculated, relative to the 1961 to 1990 mean. The anomalies from the model grid point closest to the glacier were then added to the precipitation and temperature climatologies obtained of the closest CRU CL 2. 0 grid point to force the mass balance model. All other mass balance model parameters were obtained as described in Sect. 2, and errors were propagated as described in Sect. 4.
section_name: Applied climate models

section_name: Treatment of data gaps

section_content: Since Antarctica is not covered by the CRU data sets, it is not possible to model peripheral glaciers in Antarctica using our model. Mass and surface area changes of peripheral glaciers in Antarctica were estimated by applying annually the modeled global mean specific rates of volume and surface area change.
section_name: Peripheral glaciers in Antarctica

section_content: In each region, there are glaciers that either cannot be modeled by our model because (i) there are data gaps in the digital elevation model, making the determination of either z measured terminus or z max impossible, or (ii) there are no CRU data available at the glacier's location (this is the case for some glaciers near the coast line, and on small islands). Additionally, the iterative process used to estimate the glacier's surface area at the beginning of the integration sometimes fails (see Sect. 2. 1. 3). For these glaciers, the annual, regional-mean rates of surface area and volume change were applied to estimate the region's total volume and glacier surface area change. Table 2 shows the percentages of glacier surface affected by data gaps or model failure for each region for the case of using CRU data as forcing. Numbers are very similar for the climate modelforced runs.
section_name: Other regions

section_content: Figure 17 shows the accumulated surface mass balances of the 18 regions which are directly modeled (i. e. excluding peripheral glaciers in Antarctica), converted to mm sealevel equivalent (SLE) by assuming an ocean area of 3. 62 × 10 14 m 2. All regions experienced a mass loss during the 20th century, with peripheral glaciers in Greenland being the strongest contributor to sea-level rise with almost 20 mm contribution. Remarkably, most of the mass loss here is reconstructed to have occurred during the 1930s, with an almost balanced mass budget until recently. The global, total mass loss of glaciers from 1902 to 2009 is reconstructed to be 114 ± 5 mm SLE (Fig. 18 ) 5. Rates of mass loss during the 20th century were characterized by generally faster mass loss of approximately 1. 5 mm SLE yr -1 during the first half of the century, caused by Greenland in the 1930s, Arctic Canada in the 1950s to early 1960s, and the Russian Arctic in the late 1950s and 1960s. Rates then dropped to 5 Note that the global total sums and rates shown in Figs. 18, 20, 22 and 24 include the upscaled mass balance of peripheral glaciers in Antarctica. Since these glaciers were not modeled directly, we do not have any estimate of their volume, and the right axis of these plots therefore excludes the volume of peripheral glaciers in Antarctica. Because of this, losses of more than 100 % of the indicated volume may occur. a low of around 0. 5 mm SLE yr -1 during the 1970s, and since then have been gaining speed again to currently approximately 1. 0 mm SLE yr -1. This reconstructed history is similar to the estimate of Cogley (2009), but shows higher variability, higher mass loss rates during the first half of the 20th century, and consequently a higher total reconstructed mass loss than Leclercq et al. (2011). Table 3 shows the regional, modeled surface areas and volumes in 1901 and 2009.
section_name: Results from CRU-driven model

section_content: The mass losses reconstructed using the "historical" CMIP5 data are very similar to the mass losses reconstructed from the CRU data. Even on a regional basis, the range of mass losses from the CMIP5 models captures the results from CRU data very well, with a few exceptions (see Fig. 19 ): most notably, there is only one CMIP5 ensemble member (CNRM-CM5) that produces similarly high mass losses from Greenland. Also, there is no CMIP5 ensemble member that produces as high mass losses from the low latitudes and from New Zealand as the CRU data (see also Sect. 8). Consequently, the CMIP5 ensemble mean total mass loss estimate during the CRU period is 21 mm SLE lower than the mass loss estimate from the CRU data (Fig. 20 ), even though the CRU estimate lies within the range of the CMIP5 reconstructions, both regarding the cumulative global total, as well as the rates. Table 5 gives the global total mass losses for each of the ensemble members, including the propagated model errors, during the period covered by the "historical" experiments.
section_name: Results from the model driven by data from the "historical" CMIP5 experiments

section_content: Figure 21 shows the projected regional mass losses for the RCP26, RCP45, RCP60 and RCP85 experiments until the year 2100. All regions are projected to lose glacier mass continuously throughout the 21st century, and in all regions, the ensemble range of mass losses within one RCP scenario is at least of the same order of magnitude as the difference in the means of different RCP scenarios, and in most regions, larger than the propagated model errors (not shown). This implies that climate model uncertainty, as opposed to scenario uncertainty and mass balance model uncertainty, is the most important source of uncertainty for the regional (and global) surface mass balance of glaciers in the 21st century. However, there are substantial differences between the regions in terms of uncertainty, and in terms of the fraction of glacier ice lost during the 21st century. Generally speaking, regions that contain many small glaciers (such as Central Europe with a mean glacier surface area A = 0. 31 km 2, Western Canada and US (A = 1. 03 km 2 ), Scandinavia (A = 0. 69 km 2 ), low latitudes (A = 0. 27 km 2 ) and New Zealand (A = 0. 19 km 2 )) experience higher fractional mass losses than regions that contain many big glaciers (such as Alaska (A = 3. 68 km 2 ) and Arctic Canada (A = 6. 86 km 2 )), in some cases reaching nearly complete regional deglaciation by 2100. Svalbard, even though characterized by big glaciers (A = 21. 70 km 2 ), is projected to experience a high fractional mass loss as well. However, the large mass losses are subject to very large uncertainty within each RCP ensemble. Similarly, the projections for the Russian Arctic, Greenland, and particularly Iceland are characterized by large ensemble uncertainty, probably in connection to uncertainty in the projected temperatures caused by uncertainty of future oceanic heat transport into the region. The global, total mass loss projections span the range of 82 ± 2 mm SLE to 287 ± 10 mm SLE (Fig. 22 and Table 5 ). The differences between the RCP scenarios become more obvious in the rates of mass loss towards the end of the 21st century. While the mean of the RCP26 and RCP45 ensembles indicates the onset of a stabilization by slowly declining mass loss rates from the mid 21st century onwards, the mass loss rates of the RCP60 and RCP85 projections are higher and stable approaching 2100 (Fig. 22 ). 

A few of the CMIP5 models continue the projections for the RCP26, RCP45 and RCP85 scenarios up to the year 2300 (see Table 4 ). In these long-term projections, most regions face nearly complete deglaciation in the case of the RCP85 scenarios, and most regions retain glacier ice in the cases of RCP26 and RCP45 (Fig. 23 ). The rates of mass loss approach zero towards the end of the 23rd century for all scenarios (Fig. 24 ), indicating that the glacier ice retained in the RCP26 and RCP45 is found at altitudes high enough to reach a balanced mass budget even under the increased temperatures. The projected global, total mass loss until 2300 ranges from 175 ± 2mm SLE (RCP26, MPI-ESM-LR) to 459 ± 22 mm SLE (RCP85, CSIRO-Mk3-6-0, see Table 5 ).
section_name: Results from the model driven by RCP scenario experiments

section_content: Our model neglects resolving any particular process of the surface mass and energy balance of the glaciers, in favor of applying a bulk estimate of the specific mass balance, depending on temperature and precipitation alone. But neglecting to resolve processes does not imply neglecting the impact of these processes. E. g. refreezing is an important component of the surface mass balance of some glaciers. Our model does not resolve this, but since our model is calibrated with, and validated against mass balance measurements that include effects of refreezing, the impact it has on the mass balance is included in our model. To some extend, it may be hidden in e. g. the optimal parameter values of T prec solid and T melt, but it will also be represented by the model's error. The same reasoning applies to other processes impacting the surface mass balance, such as aeolian snow transport and avalanching. While the validation demonstrates that we can be confident in the model results, therefore note, that the model is suitable for identifying the contribution of individual processes to the overall mass balance. 

While the procedure of the leave-one-glacier-out crossvalidation (Sect. 3) in principle is designed to provide an However, we prefer to assign subjective weights to different measures: because of the cumulative nature of the surface mass balance, we deem a negligible bias to be most important. Since we apply our model to projected climate change, it is then important to correctly translate changes in temperature and precipitation into mass gain or loss of the glaciers, which implies that the temperature and precipitation sensitivities of the model need to reproduce the observed variance of the surface mass balance when subjected to observed, monthly temperature and precipitation variability. Finally, we need to make sure that there is no temporal trend in the model's error, which is not measured by skill score at all. In our model validation, we completely neglect uncertainty in the measured surface mass balance values, as well as in the geodetically determined volume and surface area change measurements. This implies that the uncertainty estimates we obtain during the cross-validation, and the validation using the geodetic measurements, will rather be estimated too high than too low. 

Our model has only one glacier-specific parameter, t * that is not either given externally (such as surface area, minimum and maximum elevation, and location) or determined from climate data (such as temperature lapse rate). The obvious disadvantage of limiting the model to only one such glacier specific parameter is the relatively poor performance of the model on the glaciers on which it can be evaluated, i. e. glaciers that typically have more data available to allow for the optimization of more than one parameter. The advantage of the model design therefore only becomes apparent in the independent validation: here, the greater data-to-parameter ratio actually leads to a decrease of model uncertainty. It is therefore important to keep in mind that a model as simple as this may work comparatively well when applied to glaciers with very limited data availability, but that for glaciers for which more information is available, more complex models, potentially resolving the energy balance of the ice surface, will be more appropriate and successful. 

We find the results obtained by the model driven with CRU data in the regions low latitudes and New Zealand, i. e. the extremely high reconstructed mass losses, questionable. The validation results in Table 1 indicate that in these regions the model performance is particularly poor compared to all the other regions (but there are too few validation points available in these regions to allow for a statistically meaningful evaluation). In the low latitudes, a poor performance of the model can be expected, since the temperature index melt model applied here does not capture well the processes (e. g. sublimation) that typically are important for the surface mass balance of tropical glaciers (Mölg and Hardy, 2004; Mölg et al., 2008; Sicart et al., 2008; Winkler et al., 2009). In the case of New Zealand, the remoteness of the glaciers in terms of nearby sampled glaciers may be the cause of the weak model performance (see Fig. 12 panels c and d ). However, we find it remarkable that these two regions are the only ones where the surface mass balance of the CRU-driven model clearly lies outside the range of the model driven by "historical" CMIP5 reconstructions, particularly during the first half of the 20th century (Fig. 19 ). This may indicate that not only our model, but also potential problems in the CRU data, contribute to the weak performance. 

A somewhat surprising result from our reconstruction of the 20th century surface mass balance is that the rates of mass loss have decreased throughout most of the 20th century, after a peak around the 1930s, until very recently (Fig. 18 ). But high rates of glacier mass loss during the first half of the 20th century have been reported before: e. g. Zdanowicz et al. (2012) document high melt rates in the Canadian Arctic in the 1950s, comparable to the melt rates observed during the most recent years. While Zeeberg and Forman (2001) find a strong retreat of glaciers in the Russian Arctic occurring before 1954, they reconstruct a strongly negative surface mass balance for the Shokal'ski Glacier around 1960, coinciding with the negative surface mass balances in our reconstruction for the Russian Arctic. Finally, Chylek et al. (2006) find that the warming in Greenland between 1920 and 1930 was of similar magnitude to that during 1995 to 2005, but at a higher rate, Box et al. (2009) conclude that even the magnitude was 33 % bigger than that of the warming observed from 1994 to 2007, Fettweis et al. (2008) estimate that the surface mass loss rates of the Greenland ice sheet in the 1930s were what can be expected again only around 2100, Wake et al. (2009) conclude that, particularly in the periphery of the Greenland ice sheet, strongly negative mass balances prevailed in 1923 -1933, and Bjørk et al. (2012) ) report that many glaciers in Southern Greenland in the 1930s underwent a more rapid retreat than in recent years. Since Leclercq et al. (2011) have no length records available from the Russian and Canadian Arctic on which to base volume change reconstructions, this may also explain why, in comparison, our reconstructed past contribution of glaciers to sea-level rise is higher. It is nevertheless questionable whether the exceptionally high mass losses reconstructed from peripheral Greenland, the Canadian and Russian Arctic are real. Much of the mass loss reconstructed in these regions comes from marine-terminating glaciers, whose size is probably dynamically limited. Moreover, volume changes of floating glacier ice do not affect sea level. Since our model neither captures the ice dynamics of these glaciers, nor is able to distinguish between mass loss from ice that is afloat and land-based ice, it is very possible that the sea-level contributions from marine-terminating glaciers are overestimated in the past. In fact, Fig. 14b indicates that our model not only underestimates volume changes of calving glaciers, but that the error distribution is wider for calving glaciers than for non-calving glaciers on both sides. This implies that there are calving glaciers which have lost considerably less mass in the past than reconstructed from our model, and the lack of an adequate representation of ice dynamics in our model, as described above, may be the reason. Finally, our approach to obtain a global estimate by assuming that Antarctic peripheral glaciers behave similar to the global mean is hard to justify, and it introduces uncertainty that is hard to quantify. Upcoming studies may shed more light on the behavior of Antarctic glaciers, and may allow for better approaches to reconstructing and projecting the mass balance of the Antarctic glaciers in the future. Until then, our estimates including upscaled mass changes from Antarctic glaciers (i. e., the global sums) should be interpreted with caution.
section_name: Discussion

section_content: We have presented the construction, validation and application of a model of the global surface mass balance of glaciers from 1850 to 2300, based on observed climate data, climatemodel based reconstructions of past climate, and projections of future climate. The model is able to capture the observed surface mass balances, and the temporally integrated, geodetically measured volume and surface area changes of individual glaciers. From 1902 to 2009, based on observed climate data, glaciers are estimated to have lost 114 ± 5mm SLE of ice mass. Glaciers are reconstructed to have lost mass since the beginning of the "historical" CMIP5 experiments in 1850. Ice mass loss rates are projected to peak either around 2040 ∼ 2050 (RCP26 scenario), 2050 ∼ 2060 (RCP45 scenario), 2070 ∼ 2090 (RCP60 scenario), or around 2070 ∼ 2100 (RCP85 scenario). Independent of scenario, a new equilibrium is approached towards the end of the 23rd century, but with substantial differences between the different scenarios in the amount of ice retained -reaching from nearly complete deglaciation in some of the RCP85 scenarios (corresponding to > 450 mm SLE mass loss) to less than 50 % global ice mass loss in some of the RCP26 scenarios. Up to 2100, ensemble spread tends to be the largest source of uncertainty. Towards 2300, scenario uncertainty becomes more important. 

Supplementary material related to this article is available online at: http://www. the-cryosphere. net/6/ 1295/2012/tc-6-1295-2012-supplement. zip.
section_name: Conclusions

section_content: Acknowledgements. This work was funded by the Austrian Science Fund (FWF) : P22106-N21 and P22443-N21, and supported by the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Research Platform Scientific Computing at the University of Innsbruck. We thank Regine Hock and Graham Cogley for their comments that helped to improve the manuscript a lot. 

Edited by: V. Radic

section_name: B. Marzeion et al.: Sea-level change from glaciers

section_content: independent measure of model skill, it is not strictly independent as applied here, as we use it to determine the optimal values of the four global model parameters (Sect. 2. 2). However, the impact of optimizing four parameters using nearly 4000 mass balance measurements on the measured model skill is probably small. Most likely, it is smaller than the degrading effect the cross-validation has on measured model skill by withholding information from the model; particularly in regions with few mass balance measurements, the removal of information about one glacier (i. e. withholding the information on t * ) for the sake of evaluating the model at that glacier can be expected to have a negative impact on the measured model skill. This degradation of measured model skill may be the reason why the model performs slightly better than indicated by the propagated model error in the second, truly independent validation against geodetically measured volume and surface area changes (Sect. 5). 

In principle, it would be possible to completely automate the parameter optimization process, by maximizing the skill score of the model during the cross-validation. This would be possible since the skill score unifies different measures of model performance (i. e. correlation, bias, and variance).
section_name: B. Marzeion et al.: Sea-level change from glaciers