Global glacier volume projections under high-end climate change scenarios

. The Paris agreement aims to hold global warming to well below 2 ◦ C and to pursue efforts to limit it to 1.5 ◦ C relative to the pre-industrial period. Recent estimates based on population growth and intended carbon emissions from participant countries suggest global warming may exceed this ambitious target. Here we present glacier volume projections for the end of this century, under a range of high-end climate change scenarios, deﬁned as exceeding + 2 ◦ C global average warming relative to the pre-industrial period. Glacier volume is modelled by developing an elevation-dependent mass balance model for the Joint UK Land Environment Simulator (JULES). To do this,

precipitation and temperature lapse rates to obtain the best agreement with observed mass balance profiles. JULES is forced with an ensemble of six Coupled Model Intercomparison Project Phase 5 (CMIP5) models which were downscaled using the high resolution HadGEM3-A atmosphere only global climate model. The CMIP5 models use the RCP8.5 climate change scenario and were selected on the criteria of passing 2 o C global average warming during this century. The ensemble mean 25 volume loss at the end of the century plus/minus one standard deviation is, -64±5% for all glaciers excluding those on the peripheral of the Antarctic ice sheet. The uncertainty in the multi-model mean is rather small and caused by the sensitivity of HadGEM3-A to the boundary conditions supplied by the CMIP5 models. The regions which lose more than 75% of their initial volume by the end of the century are; Alaska, Western Canada and US, Iceland, Scandinavia, Russian Arctic, Central Europe, Caucasus, High Mountain Asia, Low Latitudes, Southern Andes and New Zealand. The ensemble mean ice loss expressed in 30 sea-level equivalent contribution is 215.2 ± 21.3 mm. The largest contributors to sea-level rise are Alaska (44.6 ± 1.1mm), Arctic Canada North and South (34.9 ± 3.0mm), Russian Arctic (33.3 ± 4.8mm), Greenland (20.1±4.4), High Mountain Asia (combined Central Asia, South Asia East and West), (18.0 ± 0.8mm), Southern Andes (14.4 ± 0.1mm) and Svalbard (17.0 ± 4.6mm). Including parametric uncertainty in the calibrated mass balance parameters, gives an upper bound global volume loss of 281.1 mm, sea-level equivalent by the end of the century. Such large ice losses will have inevitable consequences for sealevel rise and for water supply in glacier-fed river systems.

Introduction
Glaciers act as natural reservoirs by storing water in the winter and releasing it during dry periods. This is particularly vital 5 for seasonal water supply for large river systems in South Asia (Immerzeel 2013, Lutz et al. 2014, Huss and Hock 2018 and Central Asia (Sorg et al. 2012) where glacier melting contributes to streamflow and supplies fresh water to millions of people downstream. Glaciers are also major contributors to sea-level rise, despite their mass being much smaller than the Greenland and Antarctic ice sheets (Kaser et al. 2006, Meier et al. 2007, Gardner et al. 2013. Since glaciers are expected to lose mass into the twenty first century (Radic et al. 2014, Giesen and Oerlemans 2013, Slangen et al. 2014, Huss and Hock 2015, there 10 is an urgent need to understand how this will affect seasonal water supply and food security. To study this requires a fully integrated impacts model which includes the linkages and interactions between glacier mass balance, river runoff, irrigation and crop production. The Joint UK Land Environment Simulator (JULES) (Best et al. 2011) is an appropriate choice for this task because it models these processes, but is currently missing a representation of glacier ice. JULES is the land surface component of the Met Office 15 Global Climate Model (GCM), which is used for operational weather forecasting and climate modelling studies. JULES was originally developed to model vegetation dynamics, snow and soil hydrological processes within the GCM but now has a crop model to simulate crop yield for wheat, soybean, maize and rice (Osborne 2014), an irrigation demand scheme to extract water from ground and river stores and two river routing schemes; Total Runoff Integrating Pathways (Oki 1999)(TRIP) and the RFM kinematic wave model (Bell et al. 2007). The first objective of this study is to add a glacier ice scheme to JULES to 20 contribute to the larger goal of developing of a fully integrated impacts model.
The second objective is to make projections of glacier volume changes under high-end climate change scenarios, defined as exceeding 2 °C global average warming relative to the preindustrial period (2017). The Paris agreement aims to hold global warming to well below 2 o C and to pursue efforts to limit it to 1.5 o C relative to the pre-industrial period (UNFCCC 2015), however, there is some evidence that this target may be exceeded. Revised estimates of population growth suggests there is 25 only a 5% chance of staying below 2 °C and that the likely range of temperature increase will be 2.0-4.9 °C (Raftery et al. 2017). A global temperature increase of 2.6-3.1 °C has been estimated based on the intended carbon emissions submitted by the participant countries for 2020 (Rogelj et al. 2016). Therefore, in this study we make end of the century glacier volume projections, using a subset of downscaled Coupled Model Intercomparison Project phase 5 (CMIP5) models which pass 2 o C and 4 o C global average warming. The CMIP5 models use the Representative Concentration Pathways (RCP) RCP8.5 climate 30 change scenario for high greenhouse gas emissions.
The paper is organised as follows; In Section 2 we describe the glacier ice scheme implemented in JULES and the procedure for initialising the model. Section 3 describes how glacier mass balance is calibrated and validated for the present day. Section 4 presents future glacier volume projections, a comparison with other studies and a discussion on parametric uncertainty in the calibration procedure. Section 5 discusses the results, the model limitations and areas for future development. In Section 6, we summarise our findings with some concluding remarks. 5 2 Model description JULES (described in detail by Best et al. (2011)) characterises the land surface in terms of sub-grid scale tiles representing natural vegetation, crops, urban, bare soil, lakes and ice. Each grid box is comprised of fractions of these tiles with the total tile fraction summing to 1. The exception to this, is the ice tile which cannot co-exist with other surface types in a grid box.
A grid box is either completely covered in ice or not. All tiles can be assigned elevation offsets from the grid box mean which 10 is typically set to zero as default.
To simulate the mass balance of mountain glaciers more accurately we extend the tiling scheme to flexibly model the surface exchange in different elevation classes in each JULES gridbox. We have added two new surface types, glaciated and unglaciated elevated tiles to JULES (version 4.7) to describe the areal extent and variation in height of glaciers in a gridbox ( Fig. 1). Each of these new types, at each elevation, has its own bedrock sub-surface with a fixed heat capacity. These sub-15 surfaces are impervious to water, and have no carbon content, so have no interaction with the complex hydrology or vegetation found in the rest of JULES. Because glaciated and unglaciated elevated tiles have their own separate bedrock sub-surface they are not allowed to share a gridbox with any other tiles. For instance, gridboxes cannot contain partial coverage of elevated glacier ice and vegetated tiles.
JULES is modified to enable tile heights to be specified in meters above sea level, as opposed to the default option, which is 20 to specify heights as offsets from the gridbox mean. This makes it easier to input glacier hypsometry into the model and to compare the output to observations for particular elevation bands. To implement this change, the gridbox mean elevation associated with the forcing data, is read in as an additional ancillary file. Downscaling of the climate data, described in Section 2.1, is calculated using the difference between the elevation band (zband) and the gridbox mean elevation (zgbm) For the purposes of this study JULES is set up with a spatial resolution of 0.5-degree and 46 elevation bands ranging from 0 -9000m in increments of 250m. The horizontal resolution of 0.5-degree is used because it matches the forcing data used to drive the model. The vertical resolution of 250m was used based on computational cost. The vertical and horizontal resolutions of the model can be modified for any setup.
Each elevated glacier tile has a snowpack which can gain mass through accumulation and freezing of water and lose mass 30 through sublimation and melting. JULES has a full energy balance multi-level snowpack scheme which splits the snowpack into layers each having a thickness, temperature, density, grain size (used to determine albedo), and solid ice and liquid water contents. The initialisation of the snowpack properties and the distribution of the glacier tiles as a function of height is described in section 2.3. Fresh snow accumulates at the surface of the snowpack at a characteristic low density and compacts towards the bottom of the snowpack under the force of gravity. When rain falls on the snowpack, water is percolated through the layers if the pore space is sufficiently large, while any excess water contributes to the surface runoff. Liquid water below the melting temperature can refreeze. A full energy balance model is used to calculate the energy available for melting. If all the mass in 5 a layer is removed within a model timestep then removal takes place in the layer below. The temperature at each snowpack level is calculated by solving a set of tridiagonal equations for heat transfer with the surface boundary temperature set to the air temperature and the bottom boundary temperature set to the sub-surface temperature.
A snowpack may exist on both glaciated and unglaciated elevated tiles if there is accumulation of snow.
The elevation-dependent mass balance (SMBz,t) is calculated as the change in the snowpack mass (S) between successive time 10 The scheme assumes that the snowpack can grow or shrink at elevation bands depending on the mass balance, but that tile fraction (derived from the glacier area) is static with time.

Downscaling of climate forcing on elevations 15
Both glaciated and unglaciated elevated tiles are assigned heights in meters above sea level and the following adjustments are made to the surface climate in gridboxes where glaciers are present.

Air temperature and specific humidity
Temperature is adjusted for elevation using a dry and moist adiabatic lapse rate depending on the dew point temperature. First the elevated temperature follows the dry adiabat 20 where T0 is the surface temperature, γdry is the dry adiabatic temperature lapse rate ( o Cm -1 ) and ‫ݖ∆‬ is the height difference between tile elevation and the gridbox mean elevation associated with the forcing data.
If the Tz is less than the dew point temperature Tdew then the temperature adjustment follows the moist adiabat. A moist adiabatic lapse rate is calculated using the surface specific humidity from the forcing data 25 q0 is the surface specific humidity, lc is the latent heat of fusion of water at 0 o C (2.501 x 10 6 J kg -1 ), g is the acceleration due to gravity (9.8 ms -2 ), r is the gas constant for dry air (287.05 kg K -1 ) R is the ratio of molecular weights of water and dry air (0.62198) and Tv (K) is the virtual dew point temperature The height at which the air becomes saturated z is The elevated temperature following the moist adiabat is then Additionally, when Tz < Tdew, the specific humidity is adjusted for height. The adjustment is made using the elevated air temperature and surface pressure from the forcing data using a lookup table based on Goff-Gratch formula (Landolt-Bornstein 1987). The adjusted humidity is then used in the surface exchange calculation.

Longwave radiation
Downward longwave radiation is adjusted by assuming the atmosphere behaves as a black body using Stefan-Boltzmann's 10 law. The radiative air temperature at the surface Trad,0 is calculated using the downward longwave radiation provided by the forcing data LW↓z0 Where σ is the Stefan-Boltzmann constant (5.67 x 10 -8 W m -2 K -4 ). The radiative temperature at height is then adjusted Where T0 is the grid box mean temperature from the forcing data and Tz is the elevated air temperature. This is used to calculate the downward longwave radiation LW↓z at height An additional correction is made to ensure that the gridbox mean downward longwave radiation is preserved where frac is the tile fraction.

Precipitation
To account for orographic precipitation, large scale and convective rainfall and snowfall are adjusted for elevation using an annual mean precipitation gradient (%/100m) 25 where P0 is the surface precipitation, γprecip is the precipitation gradient and Z0 is the grid box mean elevation. Rainfall is also converted to snowfall when the elevated air temperature Tz is less the melting temperature (0 o C). The adjusted precipitation fields are input into the snowpack scheme and the hydrology subroutine. When calibrating the present-day mass balance, we needed to lapse rate correct the precipitation to get sufficient accumulation in the mass balance compared to observations. The consequence of this, is that the gridbox mean precipitation is no longer conserved. We tested scaling the precipitation, in a way that conserves the gridbox mean by reducing the precipitation near the surface and increasing it at height, but this did not yield enough precipitation to get a good agreement with the mass balance observations. If the model is being used to simulate 5 river discharge in glaciated catchments, then the precipitation lapse rate could be used as a parameter to calibrate the discharge.

Wind speed
A component of the energy available to melt ice, comes from the sensible heat flux which is related to the temperature difference between the surface and the elevation level and the wind speed. Glaciers often have katabatic (downslope) winds which enhance the sensible heat flux and increase melting (Oerlemans and Grisogono 2002). It is important to represent the 10 effects of katabatic winds on the mass balance when trying to model glacier melt, particularly at lower elevations where the katabatic winds speed is highest.
To explicitly model katabatic winds would require knowledge of the gridbox mean slope at elevation bands, so instead a simple scaling of the surface wind speed is used to represent katabatic winds. Over glaciated grid boxes the wind speed is where γwind is a wind speed scale factor and u0 is the surface wind speed. The simple scaling increases the wind speed relative to the surface forcing data and assumes that the scaling is constant for all heights.
Although our approach is rather crude, we found that scaling the wind speed was necessary to get reasonable values for the sensible heat flux. This is seen when we compare the modelled energy balance components to observations from the Pasterze glacier in the Alps (Greuell and Smeets 2001). The measurements consist of incoming and outgoing short and long wave 20 radiation, albedo, temperature, wind speed and roughness length at five heights between 2205m-3325m meters above sea level on the glacier. Table S6 in the Supplementary Material lists the observed and modelled energy balance components and meteorological data, for experiments with and without wind speed scaling. The comparison shows that JULES underestimates the sensible heat flux by at least one order of magnitude and the modelled wind speed is four times lower than the observations. When we increase the wind speed to match the observations there is a better agreement with the observed sensible heat flux. 25 This is because the surface exchange coefficient, which is used to calculate the sensible heat flux, is a function of the wind speed in the model.

Glacier ice albedo scheme
The existing spectral albedo scheme in JULES simulates the darkening of fresh snow as it undergoes the process of aging (Warren and Wiscombe 1980). The growth rate of the grain is an empirically derived function of the snowpack temperature.
The snow aging scheme does not reproduce the low albedo values typically observed on glacier ice, therefore a new albedo scheme is used. The new scheme is a density-dependent parameterization which was developed for the implementation in the 5 Surface Mass Balance and Related Sub-surface processes (SOMARS) model (Greuell and Konzelmann 1994). The scheme linearly scales the albedo from the value of fresh snow, to the value of ice, based on the density of the snowpack surface. The new scheme is used when the surface density of the top 10cm of the snowpack (ρsurface) is greater than the firn density (550 kgm -3 ) and the original snow aging scheme is used when (ρsurface) is less than the firn density.
10 αλ, snow is the maximum albedo of fresh snow, αλ, ice is the albedo of melting ice, ρsnow is the density of fresh snow (250 kgm -3 ) and ρice is the density of ice (917 kgm -3 ). The albedo scaling is calculated separately in two radiation bands; visible wavelengths λ = 0.3-0.7µm (VIS) and near-infrared wavelengths λ = 0.7-5.0µm (NIR). The parameters, αvis, ice , αvis snow, αnir, ice , αnir, snow, γtemp, γprecip and γwind are tuned to obtain the best agreement between simulated and observed surface mass balance profiles for the present-day (see section 3). 15

Initialisation
The model requires initial conditions for (1) the snowpack properties (2) glaciated and unglaciated elevated tile fractions within a gridbox. The location of glacier grid points, the initial tile fraction and the present-day ice mass is set using data from the Randolph Glacier Inventory Version 6 (RGI6) (RGI Consortium 2017). This dataset contains information on glacier hypsometry and is intended to capture the state of the world's glaciers at the beginning of the 21st century. A new feature of 20 the RGI6 is a 0.5-degree gridded glacier volume and area datasets, produced at 50m elevation bands. Volume was constructed for individual glaciers using an inversion technique to estimate ice thickness created using glacier outlines, a digital elevation model and a technique based on the principles of ice flow mechanics (Farinotti et al. 2009, Huss andFarinotti 2012). The area and volume of individual glaciers have been aggregated onto 0.5-degree grid boxes. We bin the 50m area and volume into elevations bands varying from 0m to 9000m in increments of 250m to match the elevation bands prescribed in JULES. 25

Initial tile fraction
The elevated glaciated fraction is where RGI_area is the area (km 2 ) at height from the RGI6, n is the tile elevation and gridbox_area (km 2 ) is the area of the gridbox. In this configuration of the model, any area that is not glaciated is set to a single unglaciated tile fraction (fracrock) with a gridbox mean elevation. It is possible to have an unglaciated tile fraction at every elevation band, but since the glaciated tile fractions does not grow or shrink, we reduce our computation cost by simply putting any unglaciated area into a single tile fraction. 5 nBands = 37 is the number of elevation bands.
We assume there is no liquid content in the snowpack by setting this to zero. The density at each level is linearly scaled with depth, between the value for fresh snow at the surface (250kgm -3 ), to the value for ice at the bottom level (917kgm -3 ).
For the future simulations the thickness and ice mass at the bottom of the snowpack comes from thickness and volume data in the RGI6. The data is based on thickness inversion calculations from Huss and Farinotti (2012) for individual glaciers which 15 are consolidated onto 0.5-degree gridboxes. The ice mass is calculated from the RGI6 volume assuming an ice density 917kgm -3. For the other layers the ice mass is calculated by multiplying the density by the layer thickness which is prescribed above.
For the calibration period, the ice mass at the start of the run (1979) is unknown. In the absence of any information about this, a constant depth of 1000m is used which is selected to ensure that the snowpack never completely depletes over the calibration period. This consists of 995m of ice at the bottom level of the snowpack and 5m of firn in the layers above. The ice content of 20 the bottom level is the depth (995m) multiplied by the density of ice.
The snow grain size used to calculate spectral albedo (see section 2.2) is linearly scaled with depth and varies between 50µm at the surface for fresh snow to 2000µm at the base for ice. The snowpack temperature profile is calculated by spinning the model up for 10 years for the calibration period and 1 year for the future simulations. The temperature at the top layer of the snowpack is set to the January mean temperature and the bottom layer and subsurface temperature is set to the annual mean 25 temperature. For the calibration period the monthly and annual temperature comes from the last year of the spin-up. Setting the snowpack temperature this way gives a profile of warming towards the bottom of the snowpack representative of geothermal warming from the underlying soil. The initial temperature of the bedrock before the spin up is set to 0 o C but this adjusts to the climate as the model spins up. We use these prescribed snowpack properties as the initial state for the calibration and future runs. 30 3 Mass balance calibration and validation

Model calibration
Elevation-dependent mass balance is calibrated for the present-day by tuning seven model parameters and comparing the output to elevation-band specific mass balance observations from the World Glacier Monitoring Service (WGMS, 2017).
Calibrating mass balance against in-situ observations is a technique which has been used by other glacier modelling studies 5 Hock 2011, Giesen andOerlemans 2013). For the calibration, annual elevation-band mass balance observations are used because there is data available for sixteen of the eighteen RGI6 regions. For validation, winter and summer elevationband mass balance is used because there is less data available.
Random parameter combinations are selected using Latin Hyper Cube Sampling (McKay, Beckman and Conover 1979) between plausible ranges which have been derived from various sources outlined below. This technique randomly selects parameter values; however, reflectance in the VIS wavelength is always higher than in the NIR. To ensure the random sampling does not select NIR albedo values that are higher or unrealistically close to the VIS albedo values, we calculate the ratio of 15 VIS to NIR albedo using values compiled by compiled by Roesch et al (2002) . The ratio VIS/NIR is calculated as 1.2 so any albedo values that exceed this ratio are excluded from the analysis. This reduces the sample size from 1000 to 198 parameter sets.
In the VIS wavelength the fresh snow albedo is tuned between 0.99 -0.7 where upper bound value comes from observations of very clean snow with little impurities in the Antarctic (Hudson et al. 2006(Hudson et al. , 1994. The lower bound represents contaminated 20 fresh snow and comes from taking approximate values from a study based on laboratory experiments of snow, with a large grain size (110 µm) containing 1680 parts per billion of black carbon (Hadley and Kirchstetter 2012). Visible snow albedos of approximately 0.7 have also been observed on glaciers with black carbon and mineral dust contaminants in the Tibetan Plateau (Zhang et al. 2017). In the NIR wavelength the fresh snow albedo is tuned between 0.85 -0.5 where the upper bound comes from spectral albedo observations made in Antarctica (Reijmer, Bintanja and Greuell 2001). We use a very low 25 minimum albedo for melting ice (0.1) the VIS and NIR wavelengths to capture dirty debris covered ice.
The temperature lapse rate is tuned between values of 4.0 -10 o C km -1 where the upper limit is determined from physically realistic bounds and lower limit is from observations based at glaciers in Alps (Singh 2001). The temperature lapse rate in JULES is constant throughout the year and assumes that temperature always decreases with height.
The wind speed scaling factor γwind is tuned within the range 1-4 to account for an increase in wind speed with height and for 30 the presence of katabatic winds. The upper bound is estimated using wind observations made along the profile of the Pasterze glacier in the Alps during a field campaign (Greuell and Smeets 2001). Table S6 in the Supplementary Material contains the wind speed observations on the Pasterze glacier. The maximum observed wind speed was 4.6 ms -1 (at 2420 meters above sea-level) while the WATCH-ERA Interim dataset (WFDEI) (Weedon et al. 2014) surface wind speed for the same time period was 1.1ms -1 indicating a scaling factor of approximately 4.
The orographic precipitation gradient γprecip is tuned between 5-25%/100m. This parameter is poorly constrained by observations therefore a large tuneable range is sampled. Tawde et al. (2016) estimated a precipitation gradient of 19%/100m for 12 glaciers in the Western Himalayas using a combination of remote sensing and in-situ meteorological observations of 5 precipitation. Observations show that the precipitation gradient can be as high as 25%/100m for glaciers in Svalbard (Bruland and Hagen 2002) while glacier-hydrological modelling studies have used much smaller values 4.3%/100m (Sorg et al. 2014) and 3%/100m (Marzeion et al (2012). The tuneable parameters and their minimum and maximum ranges are listed in Table   1.
The model is forced with daily surface pressure, air temperature, downward longwave and shortwave surface radiation, specific 10 humidity, rainfall, snowfall and wind speed from the WATCH-ERA Interim dataset (WFDEI) (Weedon et al. 2014 regions defined by the RGI6 (Fig. 2). The best regional parameter sets are identified by finding the minimum root mean square error between the modelled mass balance and the observations. 20 Figure 3 shows the modelled mass balance profiles plotted against the observations using the best parameter set for each region. The best regional parameter sets are listed in Table 2 and the root mean square error, correlation coefficient, Nash-Sutcliffe efficiency coefficient and mean bias are listed in Table 3. Nine out of the sixteen regions have a negative bias in the annual mass balance. Notably Svalbard, Southern Andes and New Zealand underestimate mass balance by 1 m.w.eq.yr -1 .
The model performs particularly poorly for the Low Latitude region which has a large RMSE (3.02 m.w.eq.yr -1 ). This region 25 contains relatively small tropical glaciers in Colombia, Peru, Ecuador, Bolivia and Kenya. Marzeion et al (2012) found a poor correlation with observations in the low latitude region when they calibrated their glacier model using CRU data. They attributed that to the fact that sublimation was not included in their model, a process which is important for the mass balance of tropical glaciers. Our mass balance model does include sublimation, so it is possible the WFDEI data over tropical glaciers is too warm. The WFDEI data is based on ERA-interim reanalysis where air temperature has been bias corrected using CRU 30 data. The CRU data comprises of temperature observations which are sparse in regions where tropical glaciers are located.
Furthermore, the quality of the WFDEI data will depend on the performance of the underlying ECMWF model. It is possible that the ECMWF model does not include glacier ice in tropical regions. The absence of ice to cool the lower atmosphere would make the grid box mean temperature too warm. In Central Europe some of the poor correlation with observations is caused by the Maladeta glacier in the Pyrenees (Fig. 4) which is a small glacier with an area of 0.52 km 2 (WGMS, 2017). When this glacier is excluded from the analysis the correlation coefficient increases from 0.26 to 0.35 and the RMSE decreases from 2.03 to 1.73 meters of water equivalent per year.

Model validation
The calibrated mass balance is validated against summer and winter elevation-band specific mass balance for each region 5 where data is available (Fig. 5). For all regions, except Scandinavia in the summer, negative Nash-Sutcliff numbers are calculated for winter and summer elevation-dependent mass balance (Table 4). The negative numbers arise because the bias in the model is larger than the variance of the observations. There are negative biases for nearly all regions implying that melting is overestimated in the summer and accumulation is underestimated in the winter.
Some, but not all, of the bias is due to the partitioning of rain and snow based on an air temperature threshold of 0 o C. The 0 o C 10 threshold is likely too low, resulting in an underestimate of snowfall. When precipitation falls as rain or snow it adds liquid water or ice to the snowpack. The specific heat capacity of the snowpack is a function of the liquid water (Wk) and ice content (Ik) in each layer (k) where Cice = 2100 JK -1 kg -1 and Cwater = 4100 JK -1 kg -1 . The liquid water content is limited by the available pore space in the 15 snowpack, therefore changes in the snowfall (ice content) control the overall heat capacity. The underestimate in the ice content reduces the heat capacity which causes more melting than observed.
Other modelling studies have used higher air temperature thresholds; 1. Increasing the temperature threshold reduces the bias slightly, therefore another explanation is that the precipitation in the WFDEI data is too low. Although we have included the variation in precipitation with height, if the gridbox mean precipitation 25 is too low, then snowfall on the elevated tiles will be underestimated. We did not bias correct the precipitation before applying the lapse rate correction like other studies have done (Marzeion et al. 2012, Huss andHock 2015). The quality of the WFDEI precipitation maybe poor because the data is constrained by rain gauge observations which are sparse in high mountains regions and often biased towards low elevation levels. Even when observations are available snowfall at higher altitudes is often difficult to accurately measure and susceptible to undercatch by 20-50% (Rasmussen et al. 2012). 30 Winter mass balance is simulated better than summer mass balance, which is seen by the lower root mean square errors for winter in Table 4. Furthermore, the biases are larger in the summer than in the winter (Table 4). It is likely that the simple albedo scheme, which relates albedo to the density of the snowpack surface, does not perform particularly well in the ablation zone. There is an improvement in the simulated summer mass balance when the glaciated area increases. This is seen by the improved correlation in Fig. 7(D) in which the validation is repeated but only grid boxes with a glaciated area greater than 500km 2 are considered. This indicates that the overestimation in melting is more pronounced for small glaciated areas than regions with a large ice extent. The overestimate in melting in Central Europe is due mainly to the Maladeta glacier as discussed in Section 3.1. This is defined as climate change that exceeds 2 o C and 4 o C global average warming, relative to the pre-industrial period (Gohar et al., 2017). Six models fitting this criterion were selected from the Coupled Model Intercomparison Project Phase 5 (CMIP5). 10 A new set of high resolution projections were generated using the HadGEM3-A Global Atmosphere (GA) 6.0 model (Walters et al. 2017). The sea surface temperature and sea-ice concentration boundary conditions for HadGEM3-A are supplied by the CMIP5 models. All models use the RCP8.5 'business as usual scenario' and cover a wide range of climate sensitivities, with some models reaching 2 o C global average warming relative to the pre-industrial period, quickly (IPSL-CM5A-LR) or slowly (GFDL-ESM2M) ( Table 5). The models also cover a range of extreme wet or dry climate conditions. This is important to 15 consider for glaciers in the central and eastern Himalaya which accumulate mass during the summer months due to monsoon precipitation (Ageta and Higuchi 1984) and because future monsoon precipitation is highly uncertain in the CMIP5 models (Chen and Zhou 2015).
The HadGEM3-A data are bias corrected using a trend preserving statistical bias method that was developed for the first Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP) (Hempel et al. 2013). This technique uses WATCH forcing data 20 (Weedon et al. 2011) to correct offsets in air pressure, temperature, longwave and shortwave downward surface radiation, rainfall, snowfall and wind speed but not specific humidity. The method adjusts the monthly mean and daily variability in the GCM variables but still preserves the long-term climate signal. The HadGEM3-A was bi-linearly interpolated from its native resolution of N216 (~60km), onto a 0.5-degree grid, to match the resolution of the WATCH forcing data which was used for the bias correction. The daily bias corrected surface fields from the HadGEM3-A are used to run JULES offline to calculate 25 future glacier volume changes. The bias correction was only applied to data up until the year 2097, which means the glaciers projections terminate at this year. A flow chart of the experimental set up is shown in Fig. 8. The HadGEM3-A climate data was generated and bias corrected for the High-End cLimate Impact and eXtremes (HELIX) project.

Regional glacier volume projections 2011-2097
Glaciated areas are divided into 18 regions defined by the RGI6 with no projections made for Antarctic glaciers because the 30 bias correction technique removes the HadGEM3-A data from this region. JULES is run for this century (2011 to 2097) using the best regional parameter sets for mass balance found by the calibration procedure (Table 2). No observations were available to determine the best parameters for Iceland and the Russian Arctic, therefore global mean parameter values are used for these regions. End of the century volume changes (in percent) are found by comparing the volume at end of the run (2097) to the initial volume calculated from the RGI6. Regional volume changes expressed in percent for low (0-2000m), medium (2250m-4000m) and high (4250m-9000m) and all elevation ranges (0-9000m) is listed in Table 6. The total volume loss over all 5 elevation ranges is also listed in mm of sea-level equivalent in Table 6 and plotted in Fig. 11. Maps of the percentage volume change at the end of the century, relative to the initial volume are contained in the Supplementary Material in Figs. S1-S7.
A substantial reduction in glacier volume is projected for all regions (Fig. 9). Global glacier volume is projected to decrease by 64±5% by end of the century, where the value corresponds to the multi-model mean ± one standard deviation. The regions which lose more than 75% of their volume by the end of the century are; Alaska (-89 ± 2%), Western Canada and US (-100 ± 10 0%), Iceland (-98 ± 3%), Scandinavia (-98 ± 3%), Russian Arctic (-79 ± 10% ), Central Europe (-99 ± 0 %), Caucasus (-100 ± 0 %), Central Asia (-80 ± 7%), South Asia West (-98 ± 1%), South Asia East (-95 ± 2 %), Low Latitudes (100 ± 0 %), Southern Andes (-98 ± 1%) and New Zealand (-88±5%). The HadGEM3-A forcing data shows these regions experience a strong warming. In most regions this is combined with a reduction in snowfall relative to the present day, which is drives the mass loss (Fig. 10). Regions most resilient to volume losses are Greenland (-31 ± 5%) and Arctic Canada North (-47 ± 3%). Both 15 these regions lose volume. In the case of Arctic Canada North snowfall increases relative to the present day which helps glaciers to retain their mass. There is a rapid loss of low latitude glaciers which has also been found by other global glaciers models (Marzeion et al. 2012, Huss andHock 2015). Our model overestimates the melting of these glaciers for the calibration period (Fig. 3)

Mass Balance Components
In this section we examine how the surface mass balance components vary with height and how this will change in the future. Figure 12 shows the accumulation, refreezing and melting contributions to mass balance averaged over low, medium and high elevations ranges for the period 1980-2000. Sublimation is excluded because its contribution to mass balance is relatively 30 small. As expected there is more melting in the lower elevation ranges and more accumulation at the higher elevation ranges.
The refreezing component, which includes refreezing of melt water and elevated adjusted rainfall, shows no clear variation with height. This is because the refreezing component can both increase and decrease with height. Refreezing can increases towards lever elevations because there is more rain and melted water. It can also decrease if the snowpack is depleted or if there is not enough pore space to hold water because previous refreezing episodes have converted the firn into solid ice. The largest accumulation rates occur in Alaska (5.3 m.w.eq.yr -1 ) and Western Canada and US (7.3 m.w.eq.yr -1 ) between 4250m-9000m and the largest melt rates are found in the Caucasus and Middle East (-7.4 m.w.eq.yr -1 ) and the Low Latitudes (-7.6 m.w.eq.yr -1 ). 5 Figure 13 shows how the global annual mass balance components vary with time for low, medium and high elevations ranges.
At the high and medium elevations accumulation, refreezing and melting decrease leading to a reduction in mass loss as glaciers disappear towards the end of the century. At high elevations mass balance is reduced from -2.2 m.w.eq.yr -1 (-177 Gtyr -1 ) during the historical period  to -0.35 m.w.eq. yr -1 (-28 Gtyr -1 ) by the end of the century (2080-2097).

Parametric uncertainty analysis
The standard deviation in the volume losses presented above are relatively small. This is because only a single GCM was used to downscale the CMIP5 data (HadGEM3-A). The uncertainty in the ensemble mean reflects the impact of the different seasurface temperature and sea-ice concentration boundary conditions, provided by the CMIP5 models, on the HadGEM3-A 15 climate. Other sources of uncertainty in the projections can arise from the calibration procedure, observational error, initial glacier volume and area, and structural uncertainty in the model physics. It is beyond the scope of this paper to investigate all the possible sources of uncertainty on the glacier volume losses. Instead we discuss the impact of parametric uncertainty in the calibration procedure in the following section.
In the calibration procedure the mass balance was tuned to obtain an optimal set of parameters for each RGI6 region, however, 20 there may be other plausible parameter sets that perform equally well (i.e. for which the RMSE between the observations and the model is small). The principle of 'equifinality', in which the end state can be reached by many potential means, is important to explore because some parameters may compensate for each other. For example, the same mass balance could be reached by increasing the wind scale factor which enhances melting or decreasing the precipitation gradient which would reduce accumulation. To identify the experiments that perform equally well, we identify where there is a step change in the gradient 25 of the RMSE for each RGI6 region. A similar approach was used by Stone et al. (2010) to explore the uncertainty in the thickness, volume and areal extent of the present-day Greenland ice sheet from an ensemble of Latin Hypercube experiments.
The step change in the RMSE is identified using the changepoint detection algorithm called findchangepts (Rebecca, Fearnhead andEckley 2012, Lavielle 2005) from the MATLAB signal processing toolbox. The algorithm is run to find where the mean of the top ten experiments (excluding the optimal experiment) changes the most significantly. For each RGI6 region 30 the step changes in the RMSE are shown in Fig. 12. JULES is re-run for each of the downscaled CMIP5 experiments and for each parameter set that is defined as performing equally well (See Table S1 in the Supplementary Material for a list of the parameters sets that perform equally well). The volume losses expressed in mm of sea-level equivalent are shown in Fig. 13. The effects of the parametric uncertainty on the volume losses varies regionally, with the largest impact found for Central Europe and Greenland. Regional volume loses including parametric uncertainty in the calibration are summarised in Table 7. Including calibration uncertainty in this way gives an upper bound of 247.3 mm sea-level equivalent volume loss by the end of the century.
Another way to explore the uncertainty in the volume projections caused the calibration procedure, is to use different 5 performance metrics to identify best parameters sets. In addition to using RMSE, we calculate best parameter sets by (1) minimising the absolute value of the bias and (2) maximizing the correlation coefficient. The best regional parameter sets are different depending on the choice of performance metric used (See Tables S2 and S3 in the Supplementary material). For twelve regions, minimising the bias results in higher precipitation lapse rates, than when RMSE values are used to select parameters. This suggests the bias in many regions is caused by underestimating the precipitation lapse rates. As discussed 10 above, this could be due to the fact the gridbox mean WFDEI precipitation was not bias corrected. Glacier volume projections are generated by repeating the simulations using these two additional performance metrics to identify best parameter sets. The uncertainty in the global volume loss when the extra performance metrics are used, is approximately double the uncertainty arising from the different climate forcings (Fig. 16, Table 7). When extra performance metrics are used, the upper bound volume loss increases to 281.1 mm sea-level equivalent by the end of the century. 15

Comparison with other studies
We compare our end of the century volume changes (excluding parametric uncertainty), to two other published studies which used the CMIP5 ensemble under the RCP8.5 climate change scenario Hock 2015, Radic et al. 2014). Other studies exist, but these include the volume losses from Antarctic glaciers which makes a direct comparison difficult (Marzeion et al. 2012, Slangen et al. 2014, Giesen and Oerlemans 2013, Hirabayashi et al. 2013). Huss and hock (2015)  Our end of the century percentage volume losses compare reasonably well to Huss and Hock (2015) for Central Europe, Caucasus, South Asia East, Scandinavia, Russian Arctic, Western Canada and US, Arctic Canada North, North Asia, Central Asia, Low Latitudes and New Zealand but are significantly higher in the Southern Andes, Alaska, Iceland and Arctic Canada 25 South. The uncertainty in our percentage volume losses are smaller than Huss and Hock (2015) because we only use a single GCM to downscale the CMIP5 experiments while Huss and Hock (2015) use 14 CMIP5 GCMs.
We estimate the end of century global sea-level contribution, excluding Antarctic glaciers, to be 215 ± 20mm which is higher than 188mm (Radic et al. 2014) and 136±23mm (Huss and Hock 2015) caused mainly by greater contributions from Alaska, Southern Andes and the Russian Arctic. These three regions are discussed in turn. 30 For the Southern Andes our estimates are approximately double (14.4mm) that of the other studies (5.8mm (Huss and Hock 2015), 8.5mm (Radic et al. 2014)). This region has the largest negative bias in the calibrated present-day mass balance (-2.87 m.w.eq.yr -1 see Table 3). To explore the effects of correcting the calibration bias on the ice volume projections, we subtract the bias values listed in Table 3 from the future annual mass balance rates. Each gridbox is assumed to have the same regional mass balance bias. The bias corrected volume losses are listed in Table S5 in the supplementary material. For the Southern Andes, the volume losses are much closer to the other studies (7.6mm) when the bias is corrected. The impact is less for the other regions where the biases are smaller. For the Russian Arctic our volume losses are higher than the other studies but that should be interpreted with caution because there were no observations available in this region to get a tuned parameter set 5 (global mean parameters where used instead). In Alaska the bias in annual mass balance is small (0.06 m.w.eq.yr -1 ) so correcting the bias has little effect on the volume loss projection for this region. Applying the bias correction increases the global volume loss from 215 ± 20mm to 223 ± 20 mm, therefore the difference between our model and the other studies cannot be explained by the bias in the calibration.
It is likely that our SLE contributions are higher than the other studies because the climate forcing data is different. The 10 HadGEM3-A model uses boundary conditions from a subset of RCP8.5 CMIP5 models with the highest warming levels.
Furthermore, the HadGEM3-A data has a higher resolution (approximately 60km) than the CMIP5 data which was used by the other two studies. This means our model should, in theory, be more accurate at reproducing regional patterns in precipitation and temperature over complex terrain which is important for calculating mass balance.
Another explanation why we predict more volume loss than Radic et al. (2014) and Huss and Hock (2015) is because there is 15 no retreat of the glacier terminus represented in the model. The only way for glaciers to reach equilibrium with climate is by melting completely. A study by Marzeion et al. (2014) showed that models predict more mass loss when the terminus elevation is fixed than when it is allowed to vary. This is because when the terminus is allowed to retreat, there will be less area available to melt. Marzeion et al (2014) found that neglecting terminus elevation changes resulted in an extra few tens of mm SLE depending on RCP scenario. Lastly, some of the differences between our study and other published projections could be due 20 models using different initial ice volumes and glaciated areas.

Discussion
The robustness of the glacier projections depends on how well the model can reproduce present-day glacier mass balance. One of the main shortcomings of the calibration and validation of mass balance is that only a single type of observations is used.
This data was used because we wanted to ensure the model could reproduce variations in accumulation and ablation with 25 height when the elevated tiling scheme was introduced. Point mass balance observations are affected by local factors such as aspect, avalanching, debris cover and there is a possibility that these local factors affect parameter sets chosen for entire RGI region. This could be improved by using observations from satellite gravimetry and altimetry, such as that described by Gardner et al (2013) to get a quantitative estimate of the model performance at the regional scales.
One of the notable differences between our study and other global glacier models is that our tuned precipitation lapse rates are 30 very high, for example, 24%/100m for South Asia West and 19%/100m for Central Asia. Other models have used lower precipitation lapse rates (1-2.5%/100m (Huss and Hock 2015), 3%/100m (Marzeion et al. 2012)), but they also bias correct precipitation by multiplying by a scale factor. This scaling factor can be considerably high. Giesen and Oerlemans (2012) found that precipitation needed to be multiplied by a factor of 2.5 to get good agreement with mass balance observations. Radic and Hock (2011) derived, through calibration of present-day mass balance, a precipitation scale factor of as high as 5.6 for Tuyuksu and Golubina glaciers in the Tien Shan. Our lapse rates are high because we do not bias correct the precipitation using a multiplication factor for the present-day. For the future GCM data the gridbox mean precipitation was bias corrected 5 using the ISI-MIP technique. The negative bias that we get when validating the present-day mass balance suggests that snowfall is underestimated in our model. A future study using this model could test whether bias correcting the precipitation before applying the lapse rate correction improves the simulated mass balance. This is the first attempt to implement a glacier scheme into JULES and so the model has many limitations. One of the key shortcomings is that glacier dynamics is not included (glacier area does not vary). The model does not simulate the retreat of 10 the glacier terminus which results in an overestimate of mass loss. Neither does the model simulate the transport of ice from higher elevations to lower elevations. These processes could be included in future work by adding a volume-area scaling scheme Bahr (1997) or a thickness parameterisation based on glacier slope (Marshall and Clarke 2000). Volume-area scaling has been used to model glacier dynamics in coarse resolution (0.5-degree) models where all glaciers in a gridbox are represented by a single ice body (Kotlarski et al. 2010, Hirabayashi, Doll andKanae 2010). The current configuration of 15 elevated glaciated and unglaciated tiles in JULES makes it well suited to a volume-area scaling model. This would be implemented by growing (shrinking) the elevated glaciated tiles if mass balance is positive (negative) at each elevation band.
In the case where the elevated ice tile grows the unglaciated tile would shrink at that elevation band or vice versa.
The volume-area scaling law has been used successfully to model the dynamics of individual glaciers (Radic et al. 2014, Giesen and Oerlemans 2013, Marzeion et al. 2012, Slangen et al. 2014) but has some limitations when applied to coarse models 20 where glaciers are consolidated into a single gridbox. The volume-area scaling law, relates volume to area using a constant scaling exponent which is typically derived from a small sample of glacier observations (Bahr et al. 1997). One of the draw backs is that the law is non-linear, meaning the exponent derived from individual glaciers would overestimate the volume of a large ice grid box such as in our model (Hirabayashi et al. 2013). Furthermore, the exponent may not accurately represent the volume-area relationship in other geographical regions. To overcome these issues a spatially variable scaling exponent 25 could be created using the newly available 0.5-degree data on volume and area contained in the RGI6.
Another shortcoming of the model is the simple treatment of katabatic winds, which is modelled by scaling the synoptic wind speed. This could be improved by parameterising katabatic winds based on the gridbox slope and the temperature difference between the glacier surface and the air temperature using the Prandtl model (Oerlemans and Grisogono 2002).
An additional drawback of the model is the coarse resolution of the gridboxes which make it unfeasible to include some process 30 which affect local mass balance such as hillside shading, avalanching, blowing snow and calving. The model could, however, be run on a finer resolution using higher resolution climate forcing data.
While this modelling projects considerable reduction in glacier mass for all mountain ranges by the end of this century, it is clear that many of the world's mountain glaciers will evolve in ways that are currently difficult to model. For instance, paraglacial processes during deglaciation lead to enhanced rock falls and debris flows from deglaciating mountain slopes and these deliver rock debris to glacier surfaces. This produces debris-covered glaciers and these are common in many mountain regions, including in Alaska, arid Andes, central Asia and in the Hindu Kush-Himalaya. Thick debris cover (decimetres to metres) limits ice ablation, (e.g., (Lambrecht et al. 2011, Pellicciotti et al. 2014, Lardeux et al. 2016, Rangecroft et al. 2016 and reverses the mass balance gradient, with comparatively higher ablation rates up glacier than at the debris-covered terminus. 5 This significantly influences glacier dynamics (Benn et al. 2005), and with inefficient sediment evacuation eventually leads to the transition from debris-covered glaciers to rock glaciers (e.g (Monnier and Kinnard 2017). In the context of continued climate warming, the transition from ice glaciers to rock glaciers may enhance the resilience of the mountain cryosphere (Bosson and Lambiel 2016). As a result, better assessment of the response of the mountain cryosphere to climate warming will depend on a clearer understanding of glacier-rock glacier relationships. 10 There are three key strengths to the JULES glacier model. Firstly, we include variations in orography within a climate gridbox which is important to calculate elevation-dependent glacier mass balance. Kotlarski et al (2010) developed a glacier scheme for the REMO regional climate model by lumping glaciers into 0-5-degree gridboxes in a similar approach to us, but they did not have a representation of subgrid orography. Instead glacier gridboxes received double the gridbox mean snowfall, glacier ice had a fixed albedo and a constant lapse rate was applied to adjust temperatures. They concluded that to reproduce mass 15 balance trends over the Alps, the scheme needed to include subgrid variability of atmospheric parameters within a gridbox. Secondly, the model uses a full energy balance scheme to calculate glacier melting. This is a more physically based approach than the widely used temperature index models, which relate melting to temperature using a degree day factor (DDF). The DDF lumps all the energy balance components into a single number meaning that the effects of changing wind speed, cloudiness and radiation on melt rates cannot be considered. Changes in solar radiation can be an important driver of melting. 20 Huss et al (2009) studied long term mass balance trends for a site in the Alps and showed that melting was stronger during the 1940's than in recent years despite more warming. This was because summer solar radiation was higher during the 1940s.
Moreover, temperature index models have been found to be less accurate with increasing temporal resolution (for example on daily time steps) (Hock 2005). Finally, the glacier scheme is coupled to a land surface model, which presents opportunities for further studies. For instance, the model could be used to investigate the impact of climate change on river discharge in 25 glaciated catchments in Asia, South America or the Arctic.

Conclusions
The first aim of this study was to add a glacier component to JULES to develop a fully integrated model, to simulate the interactions between glacier mass balance, river runoff, water abstraction by irrigation and crop production. To do this we added two new surface types to JULES; elevated glaciated and unglaciated tiles. This allows us to the calculate elevation-30 dependent mass balance which can be used to study the response of glaciers to climate change. Glacier volume was modelled by growing or shrinking the snowpack, using the elevation-dependent mass balance, but glacier dynamics was not included. Present-day mass balance was calibrated by tuning albedo, wind speed, temperature and precipitation lapse rates to obtain a set of regionally tuned parameters which are then used to model future mass balance. Winter and summer mass balances are reproduced reasonably well for regions where the glaciated area is large, however, the model performs poorly for small glaciers particularly in the summer. The fully integrated model is potentially a useful tool for the scientific community to study the impact of climate change on food and water resources. 5 The second aim of this study was to make glacier volume projections for the future under a range of high-end climate change scenarios. The ensemble mean volume loss ± one standard deviation is, -64±5% for all glaciers excluding those on the periphery of the Antarctic ice sheet. The small uncertainties in the multi-model mean are caused by the sensitivity of HadGEM3-A to the boundary conditions supplied by the CMIP5 models. Our end of the century global volume loss is 215 ± 20mm, which is higher than values reported by other studies. This is because we used a subset of CMIP5 models with the 10 highest warming levels to drive the model and glacier dynamics is not included which results in more mass loss than other studies that include dynamics. Including parametric uncertainty in the calibration procedure results in an upper bound global volume loss of 281.1 mm sea-level equivalent by the end of the century. The projected ice losses will have an impact on sealevel rise and on water availability in glacier-fed river systems.   Wind speed scale factor 1 -4 γwind -      Comparison between modelled and observed elevation-band specific mass balance for winter (grey triangles) and summer (black dots). The modelled mass balance is calculated using the tuned regional parameters from the calibration procedure. Table 4 Winter (bold) and summer (italics) number of elevation-band mass balance observations (No obs), root mean square error (RMSE), correlation coefficient (r), Nash-Sutcliffe efficiency coefficient (NS) and mean bias (BIAS).  Figure 6 Simulated and observed elevation-dependent winter mass balance when gridboxes with a glacier area of less than 100km 2 , 300km 2 and 500km 2 are excluded. The colour identifies the RGI6 regions shown in Figure 2. The RMSE, correlation coefficient and number of glaciers are listed. Figure 7 Simulated and observed elevation-dependent summer mass balance when gridboxes with a glacier area of less than 100km 2 , 300km 2 and 500km 2 are excluded. The colour identifies the RGI6 regions shown in Figure 2. The RMSE, correlation coefficient and number of glaciers are listed. Table 5 List of high-end climate change CMIP5 models that are downscaled using HadGEM3-A. The years when the CMIP5 models pass +1.5 o C, +2 o C and +4 o C global average warming relative to the pre-industrial period are shown. *No data is available for 2113 because the bias corrected data ends at 2097.   Table 6 Percentage ice volume loss, relative to the initial volume (∆V) and ice loss in mm of Sea Level Equivalent (SLE) for the end of the century (2097). Percentage volume loses are shown for low, medium and high elevation ranges as well as for all elevations. The data shows the multi-model mean ± one standard deviation. The conversion of volume to SLE assumes an ocean area of 3.618 x 10 8 km 2 . The initial area and volume from the Randolph Glacier Inventory Version 6 is listed in columns 1 and 2.  Figure 11 (A) Regional percentage volume losses at the end of the century (2097), relative to the initial volume and (B) volume losses expressed in sea-level equivalent contributions. The large red dots represent the multi-model mean and the small black dots are the individual HadGEM3-A model runs. Figure 12 Modelled annual surface mass balance components; accumulation, refreezing and melting for the period 1980-2000 for RGI6 regions. To make the figure easier to read melting is given as positive sign and sublimation is excluded because its contribution is very small. Mass balance components are averaged over low (0-2000m), medium (2250m-4000m) and high (4250-9000m) elevation ranges.
5 Figure 13 Global mass balance components for three elevation ranges. The historical period is calculated using the WDFEI data and the future period is the multi-model means of all GCMs. The bars show the averages over the time periods for accumulation, refreezing, melt and mass balance rates. Figure 14 Calibration experiments raked according the root mean square error between simulated and observed mass balance profiles for RGI6 regions. There are 198 experiments but only the top 30 have been plotted to make the figure easier to read. The red dots indicate experiments that perform equally well. Figure 15 Regional volume losses expressed in sea-level equivalent including parametric uncertainty in mass balance parameters. The solid lines show the volume loss for each downscaled CMIP5 GCM using the optimum parameter sets. The dashed lines are for runs which use other equally 'good' parameter sets based on the RMSE. Table 7 Regional ensemble mean, minimum and maximum volume losses for 2097 in sea-level equivalent (mm) when the presentday mass balance is calibrated in different ways. Columns 1-3; mass balance is calibrated by minimising the RMSE. Columns 4-6; mass balance is calibrated using an ensemble of equally plausible RMSE values. Columns 7-9; mass balance is calibrated by minimising the RMSE, minimising the bias and maximising the correlation coefficient.

Code Availability 5
The glacier scheme is included in JULES v4.7. The source code can be downloaded by accessing the Met Office Science Repository Service (MOSRS) (requires registration): https://code.metoffice.gov.uk/ The code used for this study is in https://code.metoffice.gov.uk/svn/jules/main/branches/dev/sarahshannon/vn4.7_va_scaling