Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change

Glacier mass loss is recognized as a major contributor to current sea level rise. However, large uncertainties remain in projections of glacier mass loss on global and regional scales. We present an ensemble of 288 glacier mass and area change projections for the 21st century based on 11 glacier models using up to 10 general circulation models and four Representative Concentration Pathways (RCPs) as boundary conditions. We partition the total uncertainty into the individual contributions caused by glacier models, general circulation models, RCPs, and natural variability. We find that emission scenario uncertainty is growing throughout the 21st century and is the largest source of uncertainty by 2100. The relative importance of glacier model uncertainty decreases over time, but it is the greatest source of uncertainty until the middle of this century. The projection uncertainty associated with natural variability is small on the global scale but can be large on regional scales. The projected global mass loss by 2100 relative to 2015 (79 ± 56 mm sea level equivalent for RCP2.6, 159 ± 86 mm sea level equivalent for RCP8.5) is lower than, but well within, the uncertainty range of previous projections.


Introduction
Glaciers are prominent features of many mountain and arctic landscapes. All glaciers on Earth (excluding the ice sheets, but including peripheral glaciers in Greenland and Antarctica) cover around 700,000 km 2 (RGI Consortium, 2017). Their almost universal retreat in recent decades (Zemp et al., 2019), which has been attributed predominantly to anthropogenic forcing (Marzeion, Cogley, et al., 2014;Roe et al., 2017), has made them an icon of climate change. Glaciers act as seasonal redistributors of water input into many drainage basins, modifying water availability (e.g., Bliss et al., 2014;Huss & Hock, 2018;Immerzeel et al., 2013Immerzeel et al., , 2020. Both advancing and retreating glaciers can create geohazards, such as glacier lake outburst floods, catastrophic ice collapses, and increased slope instability (e.g., Kääb et al., 2018;Richardson & Reynolds, 2000).
Glaciers store 158±41×10 3 km 3 of ice , which, if distributed equally over today's ocean area (362.5×10 6 km 2 , Cogley, 2012), and taking into account that roughly 15 % of the ice is located below present-day sea level, corresponds to 0.32 ± 0.08 m sea level equivalent (SLE). This is less than 1 % of the global ice mass of approximately 66 m SLE (Vaughan et al., 2013). Yet glaciers are contributing significantly to sea level rise on centennial time scales. A recent synthesis of available regional estimates based on glaciological, geodetic, and gravimetric data found a glacier contribution of0.77 ± 0.31 mm SLE yr −1 (0.61 ± 0.08 mm SLE yr −1 excluding the glaciers in the Antarctic and Greenland periphery) during the period 2006-2015 (Hock, Rasul, et al., 2019) to an observed total rate of global mean sea level rise of 3.16±0.37 mm yr −1 (Oppenheimer et al., 2019). Hence, the relative contribution from glaciers (excluding Greenland and Antarctic periphery) to recent rates of sea level rise is roughly 24 %.
Glaciers are expected to remain important contributors to sea level rise in the 21st century in response to anticipated global warming (Hock, Bliss, et al., 2019;Meier, 1984). Independent of the emission scenario applied, glaciers have been projected to contribute more to global mean sea level rise by the end of the 21st century relative to 1986-2005 than either the Greenland or Antarctic ice sheet. However, the uncertainty of the projected contribution of the ice sheets is larger than that of the glaciers (IPCC, 2019). Because of a disequilibrium between present-day glacier mass and concurrent climate conditions, future mass loss is unavoidable even without any further climate change (Mernild et al., 2013;Marzeion et al., 2018;Zekollari et al., 2020); however, the contribution of this disequilibrium to total projected 21st century glacier mass losses has not been quantified yet.  provided the first intercomparison of projections of global glacier mass change, based on an ensemble of 214 projections from six different glacier models, using climate projections from 25 general circulation models (GCMs) covering four RCP scenarios as boundary conditions. This intercomparison was based on previously published glacier model projections, leading to some limitations in the analysis of the data. For example, most models used the Randolph Glacier Inventory (RGI, Pfeffer et al., 2014) as the initial condition, but different versions were applied (ranging from version 1.0 to 4.0; one model used the World Glacier Inventory and upscaled this for global coverage). There were also only four GCM projections that were used by all of the glacier models. Since the projected temperature and precipitation changes vary between GCMs, the composition of the GCM ensembles of the different glacier models likely affected the projected glacier evolution. It was thus difficult to identify the origins of differences among the projections, since both initial and boundary conditions differed among the glacier models.
Like the study of , this study is part of the Glacier Model Intercomparison Project (Glacier-MIP), a Targeted Activity by the Climate and Cryosphere project of the World Climate Research Programme. The overall goal of this study is to provide updated global-scale glacier volume and area projections for 2015-2100 based on multiple glacier models, 10 GCMs, and 4 RCPs (Figure 1). For the first time, all 11 participating glacier models used prescribed initial and boundary conditions, thus forming a more coherent ensemble of projections compared to . Four of the glacier models presented here also participated in the intercomparison by . The additional seven glacier models add diversity to the ensemble by, for example, relying on flow line modeling of ice dynamics instead of scaling relations between volume and area to represent glacier geometry change, by calculating glacier melt based on energy balance instead of temperature-index melt models, and by inclusion of a gridded Earth surface model that incorporates glaciers.
In addition, we approximate the contributions to specific glacier mass balance from (i) the disequilibrium between present-day glacier geometry and present-day climate and (ii) from future warming above the present-day temperature. The ensemble of projections is also analyzed to quantify the total uncertainty and four specific contributors to the total uncertainty: differences in glacier models, differences in GCMs, differences in RCP scenarios and natural variability.

Initial and Boundary Conditions
To participate in the intercomparison, glacier modeling groups were requested to adhere to a set of prescribed initial and boundary conditions and submit data for at least one complete RGI primary region. The participating 11 glacier models are listed in Table 1. The RGI Version 6.0 was used to obtain initial conditions. For Greenland peripheral glaciers, only those with connectivity Levels 0 and 1 were included to enable a clear separation of ice sheet and peripheral glaciers (see Rastner et al., 2012, for a definition). For glacier models that require an external data set of the ice thickness distribution or volume of the glaciers, an update of the data set from Huss and Farinotti (2012) to RGI Version 6.0 was provided. For some models, it is not possible to prescribe an ice volume without creating internal inconsistencies, and initial volumes were instead computed internally. To convert ice volume to mass, a density of 900 kg m −3 was used. To convert ice mass change to SLE, we used an ocean area of 362.5 ×10 6 km 2 (Cogley, 2012).
For boundary conditions, we selected ten GCMs from the Coupled Model Intercomparison Project Phase 5 (CMIP5) ensemble (Taylor et al., 2012;Figure 1c). The choice was guided by model performance according to Walsh et al. (2018), by availability of multiple RCPs, and by the choice in the previous model intercomparison by . Participants were required to force their glacier model with projections using radiative forcing scenarios RCP2.6 and RCP8.5, and RCP4.5 and RCP6.0 if possible. These scenarios describe pathways of radiative forcing (van Vuuren et al., 2011) and are named after the radiative forcing around 2100 relative to preindustrial values in W m −2 (i.e., 2.6, 4.5, 6.0, and 8.5 W m −2 , respectively). For 7 of the 10 GCMs, projections forced by all four RCPs were available, while the other three (CanESM2, CNRM-CM5, and MPI-ESM-LR) did not have RCP6.0 projections. From each of these models' projections, we selected the r1i1p1 ensemble member. Participants were free to downscale and/or bias correct the climate data as they chose (Table S2 in the supporting information). Figure 2 shows the global mean temperature and precipitation anomalies of the GCM simulations used as boundary conditions, as well as the anomalies only over the glacier area. The latter are computed by extracting for each glacier the air temperature and precipitation of the closest grid cell and calculating the area-weighted mean. For all GCMs and RCP scenarios, glaciers experience warming that is stronger than the global mean. This amplification of warming is not just caused by the warming being stronger over land than ocean, but also because the spatial distribution of glaciers is biased toward high latitudes, exposing glaciers to Arctic Amplification (Barnes & Polvani, 2015). The temperature anomaly amplification factor experienced by glaciers (i.e., the mean anomalies over the glaciers divided by the mean anomalies globally, averaged over the period 2015 to 2100) is 1.8 for RCP2.6, 1.9 for RCP4.5, 3.0 for RCP6.0, and 1.8 for RCP8.5. The high temperature anomaly amplification factor for RCP6.0 is likely affected by the reduced GCM ensemble size for that scenario (Figure 1).     Anomalies are given globally averaged as well as averaged only over the glacier area (weighted by each glacier's area in the RGIv6). The "historical" r1i1p1 experiment of each GCM is used to extend the RCP projections into the past for calculating the reference data. Solid lines: ensemble median. Shading: 15th to 85th percentile of ensemble. A 5-yr moving average is shown for clarity.
Though not properly resolved in GCMs, glaciers may also experience amplified temperature change through elevation-dependent warming (Pepin et al., 2015;Palazzi et al., 2018). The global mean precipitation increase throughout the century is also amplified over the glaciers. However, with increasing temperatures, the solid fraction of precipitation will be reduced. On the global scale, the impact of the precipitation increase on accumulation is therefore limited (Marzeion, Jarosch, and Gregory 2014).

Glacier Models
Annual glacier area and mass change until 2100 from 11 glacier models were submitted to the intercomparison. Results were provided as regional aggregates according to the RGI primary regions (RGI Consortium, 2017). The data are available as supporting information (Marzeion et al., 2020). We refer to these models by acronyms given by the authors, or if none was given, based on the most relevant reference (following the convention of . Four models generated results for all 19 regions and three models for all regions except Antarctica and Subantarctic, while two models only included the three regions of High Mountain Asia, and two models covered only one region (Central Europe, New Zealand). Figure 1b shows the number of GCM ensemble members from each glacier model, ranging from nine (JULES) to a complete ensemble (37 members, MAR2012 and WAL2001). Table 1 summarizes the characteristics of each of the glacier models. Four of these models took part in the first GlacierMIP study (MAR2012, RAD2014, GloGEM, and WAL2001, the latter one labeled SLA2012 in , but all projections used here were updated according to the prescribed initial and boundary conditions. All of the models calculate glacier mass balance and represent glacier geometry change, but there is considerable variety in the complexity of both the mass balance and glacier geometry change parameterizations. While none of the models in Hock et al. (2019) used a prognostic ice dynamics model, there are now two models that calculate glacier geometry change using flowline models of ice dynamics (OGGM and GloGEMflow). The remainder of the models relies on various parameterizations based on scaling between area, volume, and length (Bahr et al., 1997;van de Wal & Wild, 2001) or from empirically derived relations between glacier mass and geometry change (Huss et al., 2010). Of all the models, only GloGEM takes into account frontal ablation of marine-or lake-terminating glaciers, based on the parameterization of Oerlemans and Nick (2005).
Most models compute ablation using a temperature-index melt model and accumulation using a threshold temperature to distinguish between solid and liquid precipitation (Table 1). Most models apply a bias correction on the boundary conditions obtained from the GCMs, but reference data sets and methods vary greatly. All of the models calibrate model parameters (e.g., threshold temperatures for glacier melt or solid precipitation, lapse rates and bias corrections for temperature and/or precipitation, and degree day factors), but the number and type of parameters also varies greatly between models. Typically, the model parameter sets are obtained by maximizing the agreement between model results and observations of, for example, time series of annual and seasonal area-averaged mass changes of individual glaciers from WGMS (2017), or temporally integrated multiyear glacier mass change of individual glaciers or larger glacier regions (Cogley, 2009;Dyurgerov, 2010;Gardelle et al., 2013;Gardner et al., 2013;Scherler et al., 2011;Shean et al., 2020;Zuo & Oerlemans, 1997;Zemp et al., 2019).
Almost all glacier models evaluate model uncertainty using observations not used for calibration (WAL2001 estimates uncertainties based on parameter sensitivities, see Slangen & van de Wal, 2011). However, since every model uses different types of data, periods, and measures to evaluate model performance, a direct comparison of model quality is hampered. Consequently, there is no possibility to rank glacier models by their performance, or to introduce a weighting scheme when calculating ensemble mean values. Table S2 summarizes the calibration and evaluation procedures of the glacier models. We refer to the original publications (listed in Table 1) for details.
When specifying results in SLE, we assume that all glacier mass change contributes to sea level change, ignoring potential effects such as replenishing of groundwater, storage of melt water in lakes, or ocean area change from isostatic effects and shoreline migration. Huss and Hock (2015) found a reduction of 11 % to 14 % of the sea level contribution until 2100 from glaciers by accounting for ice grounded below sea level, which does not contribute to sea level rise when melted. Farinotti et al. (2019) estimate that roughly 15 % of the global glacier ice mass is below present-day sea level. GloGEM and OGGM are able to account for this, but for better comparability, and since there are interactions between the treatment of frontal ablation and the estimate of glacier volume (Recinos et al., 2019) that impede the transferability of the estimate of Huss and Hock (2015) to other models, we neglect this effect here.

Glacier Mass and Area Projections
Following , all results are shown for the period 2015 to 2100, although all glacier models started their simulations prior to 2015 (Table 1). The global ice mass estimate in 2015 varies between glacier models that cover all regions, ranging from 391 mm SLE in GLIMB to 551 mm SLE in RAD2014. While initial glacier volumes were provided to all modeling groups, differences in each model's starting year of the glacier projections (Table 1) cause the differences in ice volume since ice mass evolves uniquely in each model prior to 2015. In addition, some models estimated initial volumes internally. The mean of global ice mass estimates of the glacier model ensemble in 2015 (439 mm SLE) is higher than a consensus estimate of 382 ± 99 mm SLE obtained combining different methods of ice thickness estimation , but within the uncertainty range. Figure 3 shows the modeled specific mass balance rates (i.e., rates of mass change per unit area) for all regions and all four RCPs from 2015 to 2100. Under RCP2.6, most regions either exhibit relatively constant (but negative) or increasing mass balances, with some regions approaching an equilibrium by the end of the 21st century, albeit with greatly reduced glacier area and volume (Central Europe, Low Latitudes, and New Zealand). Under RCP8.5, all regions are projected to have increasingly negative mass balance rates, which is especially true for regions affected by Arctic Amplification (e.g., Svalbard and Russian Arctic), where the specific mass balance rate is projected to drop below −4000 kg·m −2 ·yr −1 by the end of the 21st century.

Evolution of Glacier Area
Independent of scenario and region, projections show decreasing glacier areas over the 21st century in response to the negative mass balances (Figure 4). In regions that are characterized by relatively small glaciers (e.g., Western Canada and US, Scandinavia, North Asia, Central Europe, Low Latitudes, Caucasus, and New Zealand), relative area loss is greater than in regions dominated by large glaciers. While glacier geometry changes are delayed with respect to glacier mass balance, in some of these regions glacier area is projected to stabilize toward the end of the 21st century under RCP2.6, as glaciers approach an equilibrium. On the other hand, some combinations of glacier model and GCM project almost complete glacier loss in various regions under scenario RCP8.5 ( Figure S8). On the global scale, most projections estimate more than half of the glacier area remains by 2100.

Mass Change Rates
Projected mass change rates depend strongly on the emission scenario and region ( Figure 5). The stronger the radiative forcing, the higher and later the peak mass loss rates. On the global scale, under RCP2.6 mass change rates peak in the first half of the 21st century slightly above 1 mm SLE yr −1 (ensemble mean) before slowly declining to roughly 0.7 mm SLE yr −1 by 2100. Under RCP8.5, the maximum occurs around the end of the 21st century, with ensemble mean rates well exceeding 2 mm SLE yr −1 , and some ensemble members exceeding 3 mm SLE yr −1 . The peak mass loss rate is caused by a combination of stabilizing specific mass balance rates (especially relevant under RCP2.6, Figure 3) and decreasing glacier volumes (particularly relevant under RCP8.5). Depending on the initial regional ice mass, different regions therefore show distinctly different evolutions of mass change rates. Regions containing many large glaciers (Antarctic and Subantarctic, and Arctic Canada North) show increasing mass loss rates throughout the 21st century under RCP8.5 (i.e., the peak can be expected later than 2100) since the highly negative specific mass balances ( Figure 3) apply over relatively large remaining glacier areas ( Figure 4). Regions characterized by small glaciers (e.g., Western Canada and US, Central Europe, Caucasus, and Low Latitudes) show declining mass loss rates throughout the 21st century, independent of RCP scenario, since either the remaining ice area quickly declines (dominant in high-and medium-emission scenarios) or the glaciers stabilize as they retreat to higher elevations (dominant for RCP2.6). The other regions show patterns that are qualitatively similar to the global changes with initially increasing rates followed by decreasing rates but peak mass loss rates, and their timing, differ among regions.

Evolution of Glacier Mass
The evolution of glacier mass relative to 2015 is shown in Figure 6 and the corresponding cumulative contribution to sea level change in Figure 7. Almost all combinations of glacier model, RCP and GCM, show substantial mass losses in every region by the end of the century (Figure 7). The single exception is JULES, which projects stable or even increasing glacier mass in some regions under RCP2.6 and RCP4.5 ( Figures  S17 and S18). The ensemble median global mass loss (±90 % confidence interval) by 2100 relative to 2015 ranges from 18 ± 13 % under RCP2.6 to 36 ± 20 % under RCP8.5. This implies that even under scenario RCP8.5, roughly two thirds of the global glacier mass remains in 2100 (Figure 6), indicating the strong potential of glaciers to continue to contribute to sea level rise in the 22nd century. By 2100, the projected ensemble median global mass losses correspond to sea level contributions of between 79 ± 56 mm SLE (RCP2.6) and 159 ± 86 mm SLE (RCP8.5).
The global mass loss is dominated by the contributions from Alaska, the Canadian Arctic, Greenland Periphery, Svalbard, the Russian Arctic, and Antarctic and Subantarctic (all having ensemble mean contributions greater than 10 mm SLE between 2015 and 2100 under RCP8.5, Figure 7) due to their considerable initial of all regional ensemble medians, weighted by glacier area of the respective region. Regional panels are sorted by descending regional glacier mass here and in all following plots. Size of ensemble and number of glacier models (in brackets) are shown in the lower left corner. Figures S1-S4 show the results for the individual glacier models.    glacier mass. The regions that contain only a small amount of the initial global glacier mass, such as Western Canada and US, Central Europe, Caucasus, and Low Latitude, have the most relative mass loss with ensemble mean mass losses greater than 85 % between 2015 and 2100 under RCP8.5 ( Figure 6). However, their total mass loss and corresponding sea level contribution is small (Figure 7) due to their small glacier volume.
Remarkably, the projected global mass loss until 2100 is very similar for RCP6.0 and RCP4.5, even though global mean temperature increase is greater for RCP6.0 than for RCP4.5 at the end of the 21st century ( Figure 2). Since the ensemble for RCP6.0 is significantly smaller than the ensemble for RCP4.5, a direct comparison is difficult. However, up until roughly 2070, temperature change at the glacier locations is higher in RCP4.5 than in RCP6.0, following radiative forcing, which is higher in RCP4.5 than RCP6.0 until middle of the century (van Vuuren et al., 2011). While the difference in mass loss rates between RCP4.5 and RCP6.0 ( Figure 5) is qualitatively similar to the difference in radiative forcing (higher mass loss rates for RCP6.0 from the middle of the 21st century onward), the mass loss initially is higher under RCP4.5. The accumulation of the mass losses then results in a dependency of glacier mass on the emission pathway, as opposed to the instantaneous greenhouse gas concentration.

Contribution of Glacier Disequilibrium to Mass Change Evolution
The annual specific mass balance of a glacier that is not in steady state can be interpreted as combining two signals: the instantaneous mass balance response to concurrent atmospheric forcing if the glacier was in steady state, and the disequilibrium between glacier geometry and longer-term atmospheric conditions. The latter occurs because the glacier adjusts its geometry (including retreat or advance) only slowly (often over decades or longer) in response to a change in glacier mass, which in contrast occurs instantaneously to changing atmospheric forcing. Hence, a glacier subjected to a step-like warming will be larger and extend to lower elevations than the new steady-state geometry that the glacier eventually will attain.
To approximately disentangle the two signals, we calculate a linear regression of the form where ΔM is the annual specific mass balance in kg m −2 and T is the annual temperature anomaly relative to the mean of 1995 to 2014 in K, computed as area-weighted temperature anomaly of each glacier's closest grid cell (Figure 2). The "historical" r1i1p1 experiment of each GCM is used to extend the RCP projections into the past. The is a measure of the temperature sensitivity of the mean specific mass balance, M diseq is a measure of the disequilibrium between the glacier and the atmospheric conditions corresponding to the mean of 1995 to 2014, and is an error term. In this analysis, we neglect any effects from precipitation anomalies.
First, for every glacier model, GCM simulation and RCP, we extract the time series of ΔM and T averaged over all RGI regions covered by that glacier model. Thus, the domain which the annual averages of these variables refer to ranges from one single RGI region (GloGEMflow, AND2012) to all 19 RGI regions (Table 1 indicates which glacier models cover which regions). For each glacier model and each year, we then perform a linear regression (equation (1)) between all values of ΔM on ΔT, regardless GCM and RCP. Hence, we generate a time series of annual and ΔM diseq for each glacier model. We combine all GCM ensemble members and all RCPs of the given year to ensure a sufficiently large sample for the regression. For comparison, we also performed the regression for each RCP separately. Results were similar but noisy due to smaller sample size.
For each region and year we compute the correlation and significance and discard all years with (p > 0.1). Particularly, early in the simulation period (i.e., close to 2015), the correlations are sometimes insignificant, since temperature anomalies are small while we ignore precipitation changes. Figure 8 shows the area-averaged M diseq and of for each glacier model as a function of time. The multimodel median of ΔM diseq , averaged over 2005 to 2014, is −533 kg m −2 yr −1 and the standard deviation between the glacier models is 585 kg m −2 yr −1 .As the glaciers lose mass and retreat in response to increasing temperature, the disequilibrium term approaches zero in the second half of the 21st century. This corresponds to a state where glaciers would be in equilibrium if they experienced 1995 to 2014 mean temperatures. In some models ΔM diseq becomes slightly positive toward the end of the 21st century, which corresponds to a state where glaciers would gain mass if they experienced 1995 to 2014 mean temperatures. The only model that does  (1)). Note that models representing only a few regions (Table 1) will show greater variability than global models. Annual values where no significant correlation between ΔM and ΔT was found (p>0.1) were discarded.
A 5-yr moving average is shown for clarity.
not show a clear tendency to equilibrate is WAL2001, which does not account for changes in glacier terminus elevation in response to glacier mass change and can reach equilibrium only through complete loss of glaciers. The temperature sensitivity is relatively constant in time for most models (Figure 8b) but differs strongly between the individual glacier models, with temporal means between−164 kg·m −2 ·yr −1 ·K −1 (JULES) and−1379 kg·m −2 ·yr −1 ·K −1 (AND2012, likely affected by representing just New Zealand, see Table 1). Of the global models, RAD2014 and GloGEM have the most negative values of −516 and−445 kg·m −2 ·yr −1 ·K −1 , respectively. The temporal mean of the multimodel median is −329 kg·m −2 ·yr −1 ·K −1 , the standard deviation of the glacier models' temporal means is 334 kg·m −2 ·yr −1 ·K −1 .

Method
Uncertainty in projected glacier change originates from five different sources: (1) glacier model uncertainty, including uncertainty from any downscaling of atmospheric conditions internal to the glacier model, which causes any two glacier models to project different glacier evolution even if the boundary and initial conditions are identical; (2) climate model uncertainty, which causes two GCMs to respond differently to identical radiative forcing, and which enters the glacier model projections through the boundary conditions (when calculating the surface mass balance); (3) scenario uncertainty, which reflects the uncertainty of the future radiative forcing affecting the GCM projections; (4) internal climate variability, that is, natural fluctuations of climate that arise without any changes in the radiative forcing of the climate system; and (5) uncertainties in the glacier inventory, such as initial glacier volume and area.
Scenario uncertainty (3) is conceptually different from the other sources of uncertainty, since instead of a lack of knowledge about it, or approximations to the natural system, it is associated with future societal decisions. However, it contributes to uncertainties in projecting glacier mass just as any other sources of uncertainty, and is thus treated here in the same way. We neglect uncertainties in the glacier inventory (5), since all models used the RGIv6.0, and since no alternative, globally complete inventories exist that would allow a sampling of this uncertainty. We assume the remaining four uncertainties to be independent, neglecting any interaction terms in the equation for the total variance across the ensemble 2 tot . Thus, 2 tot is given by where 2 gla is the variance across different glacier models, 2 GCM is the variance across different GCMs, 2 RCP is the variance across different RCPs, and 2 nat is the variance caused by natural variability. In our ensemble, we sample the first three of these uncertainties well through the inclusion of eleven glacier models, ten GCMs and four RCPs. However, our ensemble does not include a sample where the natural 10.1029/2019EF001470 variability is isolated; that is, we do not have multiple realizations where the combination of GCM, glacier model and RCP is invariant. For each GCM only the r1i1p1 realization was used in the glacier projections. Therefore, we adopt an alternative approach: The differences in total radiative forcing of the RCP scenarios are very small in the first few years. For example, in 2015, assuming a natural forcing of 280 ppm CO 2 equivalent, the difference in total anthropogenic forcing between RCP2.6 (127.97 ppm CO 2 equivalent) and RCP8.5 (128.90 ppm CO 2 equivalent) is 0.7 % (Meinshausen et al., 2011). It is thus a reasonable approximation to completely assign the differences across RCPs in the first few years of the projections to natural variability, as long as glacier model and GCM are invariant. This is equivalent to the assumption that differences between two RCP runs of one GCM are caused only by different phases of natural variability in the first couple of years of the projection. Based on this approximation, we estimate 2 nat by 2006, … , 2015) ( 3) where i = 1, … , N gla−GCM indicates the available combinations of glacier model and GCM for which the variance across RCPs is calculated. A temporal mean is calculated over the years 2006 to 2015 and assumed to be representative for the whole period 2015 to 2100. We are thus not able to estimate a temporal change of natural variability. Since the majority of glacier models do not model every region, N gla−GCM varies between regions and ranges from 37 (Antarctic and Subantarctic) to 78 (the three High Mountain Asia regions), compared to a total number of 110 available combinations of glacier model and GCM.
We estimate 2 RCP in every year t by 2 taking into account that 2 RCP,i (t) is influenced not only by the choice of RCP but also by natural variability. We estimate 2 GCM by 2 analogous to 2 RCP . N gla−RCP varies between 11 (Antarctic and Subantarctic) and 30 (the three High Mountain Asia regions) across regions, compared to a total of 44 available combinations of glacier model and RCP.
Finally, 2 gla is approximated by Since the glacier models are forced with identical natural variability in each of the N GCM−RCP subensembles, 2 nat does not have to be corrected for. N GCM−RCP is 37 for all regions (three GCMs do not provide projections for RCP6.0).
Our approach to partition the total ensemble uncertainty is similar to the one used by Hawkins and Sutton (2009) for the CMIP3 ensemble, except that (i) we have one additional source of uncertainty being the glacier models itself; (ii) we estimate natural variability based on the assumption that small differences in anthropogenic forcing lead to small changes in climate, instead of fitting polynomials to separate systematic climate change from natural variability; and (iii) we do not assign weights to ensemble members based on past performance. However, our partitioning gives the opportunity to test whether the approximations mentioned above, the incomplete sampling of the ensemble, and the assumption of negligible interactions between glacier model, GCM, and scenario, have an impact on the partitioning of total uncertainty. For this purpose, we compare tot to the ensemble standard deviation ensemble obtained without any deconstruction. Figure 9 shows the temporal evolution of total uncertainty 2 tot and its four components of specific mass balance rates averaged over every RGI region and the globe. In all regions, total uncertainty and its individual 10.1029/2019EF001470 components generally increase throughout the 21st century. At the regional scale, the magnitude of projection uncertainty associated with natural variability is between 170 and 500 kg m −2 yr −1 (Figure 9). The larger the RGI region's glacier extent, the smaller the influence of natural variability tends to be. On the global scale, natural variability contributes only about 100 kg·m −2 ·yr −1 . This indicates that a spatial decorrelation of the natural variability leads to compensating effects at larger scales. Uncertainties associated with glacier models and GCMs are of similar order of magnitude in many regions and larger than the uncertainty from natural variability, and both tend to grow slightly over time. However, in some regions (most notably, the Low Latitudes and Antarctic and Subantarctic), the glacier model uncertainty initially is particularly large and then decreases. These are the regions where the glacier model disagreement on M diseq (equation (1)) is particularly large. Disagreement on the temperature sensitivity of glaciers, and on the response of glacier geometry to mass change, is then contributing to the growth of the glacier model uncertainty over the 21st century. The uncertainty associated with the RCP scenario is by construction (following equation (4)) zero during 2005 to 2014. It then grows almost linearly over time, to values of roughly between 750 and 2000 kg·m −2 ·yr −1 at regional scales, and approximately 600 kg·m −2 ·yr −1 at the global scale at the end of the 21st century. Figure 10 shows the relative contribution of each source of uncertainty to total uncertainty. Both at regional scales and the global scale, emission scenario uncertainty becomes the dominant source of uncertainty of projected mass balance in the second half of the 21st century. Only in Iceland, a region where GCMs projections disagree relatively strongly due to uncertainty in the evolution of North Atlantic oceanic heat transport (Collins et al., 2013), uncertainty from the GCMs remains the dominant cause of total uncertainty throughout almost the entire 21st century. In relative terms, both glacier model uncertainty and the impact from natural variability diminish over time ( Figure 10). The latter tends to decrease at a lower rate in some of the smaller regions (e.g., Caucasus and New Zealand).

Uncertainty of the Specific Mass Balance
We compute the total ensemble uncertainty from the sum of the components ( tot , equation (2)) and also directly from the entire ensemble ( ensemble ). ensemble ≠ tot indicates interactions between sources of uncertainty. On the global scale, the similarity of both quantities corroborates our assumption that interactions between the different sources of uncertainty can be neglected (Figure 9). On the regional scale, our method tends to overestimate the ensemble uncertainty. Since the glacier model ensembles are different between global and regional scales, a direct comparison is difficult.

Uncertainty of the Evolution of Glacier Mass
In contrast to the general increases in total uncertainty of specific mass balance with time ( Figure 10), there are considerable regional-scale variations in the evolution of total uncertainty in the rates of mass loss ( Figures S21 and S22). Total uncertainty decreases over time in those regions where, independent of the scenarios, glaciers will strongly shrink throughout the 21st century, such that the remaining glacier area and mass become small at the end of the 21st century. However, uncertainty of rates of mass loss increases in those regions where glaciers are particularly large (e.g., Arctic Canada North, Russian Arctic, and Antarctic and Subantarctic), and this growth of uncertainty is driven mostly by emission scenario uncertainty, and to a lesser degree by GCM uncertainty. The remaining regions show an initial growth of uncertainty, caused mostly by increasing GCM and RCP uncertainty, and later on a decline in uncertainty caused by the shrinkage of remaining ice. On the global scale, for mass change rates, increasing uncertainty from the emission scenario is almost balanced by decreasing uncertainty from the glacier models and natural variability, in both cases reflecting the reduced glacier area and mass. Total uncertainty of mass loss rates therefore remains almost constant, showing only a slight increase from roughly 1 mm SLE yr −1 at the beginning of the 21st century to 1.2 mm SLE yr −1 at the end of the century.
For temporally accumulated glacier mass change over 2015 to 2100, the global total uncertainty grows almost linearly, reaching approximately 70 mm SLE in 2100 ( Figure S23). This growth is driven by accumulating uncertainty from glacier models and GCMs, and from roughly the mid 21st century onward by the scenario uncertainty that is accelerating. Also, for most individual regions, absolute uncertainty of accumulated mass change is growing over time, except for those regions that will experience a nearly complete loss of glaciers within the 21st century (most notably, Western Canada and US, Central Europe, Caucasus, and Low Latitudes).
In terms of relative contributions to total uncertainty (Figure 11), there are considerable regional differences: while scenario uncertainty is playing a growing role everywhere typically from midcentury onward, some  (6)) to the total uncertainty ( tot , equation (2)) of projected specific mass balance rate in kg·m −2 ·yr −1 . Also shown is the uncertainty obtained directly from the entire ensemble without any decomposition ( ensemble ). regions are dominated by glacier model uncertainty throughout the century (e.g., Western Canada and US, North Asia, and Low Latitudes), others by GCM uncertainty (particularly Iceland). While natural variability plays an important role in small regions that experience little compensation associated with spatial decorrelation of variability within the region (e.g., New Zealand, Caucasus, and Central Europe), on the global scale its impact is almost negligible toward the end of the 21st century. On the global scale, at the end of the 21st century, scenario uncertainty is the dominant cause of total uncertainty, roughly as important as the joint uncertainty from glacier models and GCMs.

Discussion
The ensemble median projections of global mass loss presented here are slightly lower than those in  but well within their uncertainty range (Table 2). In contrast to  who computed global-scale means and uncertainties only from the glacier models that covered all 19 RGI regions, here we compute the global uncertainties based on all available regional projections, thus also considering the glacier models that only compute a few regions. We use the standard deviation of the regional ensemble, i , where i is an RGI region, as our measure of uncertainty of the regional projections. When calculating the uncertainty of the global values from regional values, the covariance in the regional uncertainties needs to be taken into account. Hence, we compute the uncertainty of the global values global as where i, j = 1, 2, … , 19 are the RGI regions, i (t) is the uncertainty in region i, and cov i,j (t) is the covariance of uncertainties in regions i and j.
To evaluate the robustness of this approach, we compare the resulting total uncertainties for the four glacier models that compute all 19 RGI regions to those obtained if the uncertainties were computed directly from the globally averaged projections. Uncertainties obtained by the two different methods agree well (not shown). To obtain the 90 % confidence interval, we multiply the obtained standard errors with 1.645 (Table 2).
For direct comparison with Hock et al. (2019), we recalculated their results using the approach adopted here. For glaciers outside of Greenland and Antarctica, differences in relative mass and area loss are negligible except for the RCP6.0 scenario. This indicates that these differences can be attributed to the small ensemble size both for scenario RCP6.0 and for the Antarctic and Subantarctic region. Additionally, mass balance observations in the Antarctic and Subantarctic region are not dense enough to constrain glacier models well, which leads to great glacier model uncertainty ( Figure 11). The uncertainties are higher in this study than in the study by , reflecting the greater diversity of glacier models.
The reduced mass loss projections can to some degree be attributed to the addition of new glacier models to the ensemble. Compared to , three glacier models with global (GLIMB) and nearly global (JULES and OGGM) coverage were added to the GlacierMIP ensemble. Two of these models (GLIMB and JULES) project mass balances generally less negative than the ensemble mean ( Figures S1 to S4). For GLIMB, the temperature sensitivity is close to that of the ensemble median (Figure 8b), but the glaciers are closer to equilibrium with climate than the ensemble median ( Figure 8b). For JULES, the disequilibrium is close to the ensemble median, but the temperature sensitivity is considerably lower (see below for a further discussion). The impact of additional glacier models for individual regions is limited on the global scale but can be expected to be greater on the regional scale. To benefit from the diversity of glacier models in the overall ensemble, global values should therefore be derived as the sum of regional values, instead of only taking truly global models into account.
In Figure 3, we compare the global mean specific mass balance obtained from the four glacier models covering all 19 RGI regions to the average of all regional ensemble means based on all 11 models, weighted by glacier area of the respective region. This comparison shows that the mean of the four models with global coverage (Table 1) project a more negative mean mass balance than the mean of the entire ensemble, which is dominated by the seven glacier models with almost global coverage. This indicates that increasing the size of the glacier model ensemble can have a noticeable effect on the results. Therefore, even though the number of global glacier models has increased since the publication of the first intercomparison , development of new global-scale glacier models will result in added value.
The difference in how glacier models translate climate change into glacier mass change, as evident from Figure 8b, is one of the major sources of the uncertainty in projections. However, part of the model differences is also associated with the bias correction and/or downscaling scheme each model applies to obtain more localized climate conditions than are provided by the GCMs. Almost all of the glacier models use some form of bias correction. Without bias correction, positive temperature biases would erroneously erase glaciers in some locations. At the same time, negative temperature biases would not erroneously create glaciers in new locations, since no initial conditions for these glaciers exist (the same argument holds for precipitation biases). This means, even if a GCM had no global mean temperature or precipitation bias, glacier models would systematically project too little global ice mass without local bias correction.
Downscaling schemes are typically based on lapse rates for temperature. Precipitation downscaling often includes scaling factors and sometimes a lapse rate. Downscaling mimics small-scale atmospheric processes that are relevant for the glacier mass balance, but not appropriately resolved by the GCMs. These processes include orographic precipitation, windblown snow and avalanches, and small-scale feedbacks between the ice surface and air flow. These downscaling schemes are often intrinsically linked to the glacier model calibration, and their effect is thus hard to separate from differences in glacier model physics. For example, the glacier models' sensitivities of the mass balance to temperature anomalies shown in Figure 8b to some degree reflect how the models treat precipitation. If a model's downscaling increases precipitation by including a precipitation correction factor>1, which is often necessary to achieve agreement between modeled and observed accumulation rates (e.g., Huss & Hock, 2015;Giesen & Oerlemans, 2012;Marzeion et al., 2012;Rounce et al., 2020), calibration with observed mass balances will lead to a greater temperature sensitivity than if there was no correction of precipitation. Thus, prescribing common downscaling routines would necessitate a recalibration of the glacier models (potentially deteriorating their performance), creating a new obstacle to comparability.
It is therefore difficult to further partition the uncertainty of projections associated with the glacier model ensemble (Figures 1-11) into uncertainties associated with downscaling, mass balance calculations, and representation of glacier geometry change in a quantitative way. It is possible, however, to qualitatively address the sources of disagreement between glacier models. Apart from the downscaling approach, differences in the approach to the mass balance calculation are responsible for the diversity of sensitivities shown in Figure 8b. For example, the low sensitivity of JULES is best understood as a consequence of its energy Earth's Future 10.1029/2019EF001470 balance approach to the mass balance calculation: relatively small changes to both downwelling long-and shortwave radiation in some regions (Shannon et al., 2019) lead to less negative projected mass balances from JULES compared to more negative mass balances from temperature-index melt models.

Conclusions
The greatest source of uncertainty in projected glacier mass at the end of the 21st century is the lack of knowledge on future emissions, similar to global mean temperature change (Collins et al., 2013). Emission scenario uncertainty is followed by uncertainty of GCMs and glacier models. On the global scale, the impact of natural variability is very small compared to the other sources of uncertainty.
While uncertainty related to the glacier model ensemble initially is dominant and slightly grows during the 21st century (Figure 9), the uncertainties associated with RCP scenario and GCM grow faster, such that the relative importance of glacier model uncertainty shrinks over time ( Figure 11). This implies that the longer the time scale considered, the smaller is the potential to narrow uncertainty by improving glacier models. However, reducing glacier model uncertainty can have strong benefits for reducing overall uncertainty for projections within the first half of the 21st century. For such projections with a lead time of a few decades, the glacier models' estimated imbalance between glacier geometry and concurrent climate ( Figure 8) needs to be better constrained. We therefore recommend further model intercomparisons that extend to glacier reconstructions, in order to facilitate standardized tests of model performance in comparison to observed glacier changes. Such an intercomparison would also allow including estimates of committed, but not yet realized mass change, in order to investigate the imbalance between glacier geometry and concurrent climate (Mernild et al., 2013;Marzeion et al., 2018) in greater detail.
The projected glacier mass loss during the 21st century is lower than in a previous intercomparison . This reduction of projected mass loss is caused by the inclusion of two new glacier models that have a lower overall climate sensitivity than the previous glacier model ensemble. However, the difference is generally small and well within the uncertainty range of projections, and the projected contribution of glaciers to sea level rise remains comparable to that of the ice sheets in Greenland and Antarctica presented in Oppenheimer et al. (2019).