2.1 Downscaling of climate forcing on elevations

Both glaciated and unglaciated elevated tiles are assigned heights in metres above sea level and the following adjustments are made to the surface climate in grid boxes in which glaciers are present.

2.1.1 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:

$$\begin{array}{}\text{(3)}& T_z=T_{\mathrm{0}}-{\mathit{\gamma }}_{\mathrm{dry}}\mathrm{\Delta }z,\end{array}$$

where T₀ is the surface temperature, γ_dry is the dry adiabatic temperature lapse rate (^∘C m⁻¹), and Δz is the height difference between tile elevation and the grid box mean elevation associated with the forcing data.

If T_z is less than the dew point temperature T_dew then the temperature adjustment follows the moist adiabat. A moist adiabatic lapse rate is calculated using the surface specific humidity from the forcing data.

$$\begin{array}{}\text{(4)}& {\mathit{\gamma }}_{\mathrm{moist}}={\displaystyle \frac{\left(\frac{g(\mathrm{1}+lc\cdot q_{\mathrm{0}})}{r\cdot T_{\mathrm{v}}(\mathrm{1}-q_{\mathrm{0}})}\right)}{\left(\frac{C_p+lc\cdot \mathrm{2}\cdot q_{\mathrm{0}}\cdot R}{r\cdot T_{\mathrm{v}}\mathrm{2}(\mathrm{1}-q_{\mathrm{0}})}\right)}}\end{array}$$

q₀ is the surface specific humidity, l_c is the latent heat of fusion of water at 0 ^∘C (2.501×10⁶ J kg⁻¹), g is the acceleration due to gravity (9.8 m s⁻²), r is the gas constant for dry air (287.05 kg K⁻¹), R is the ratio of molecular weights of water and dry air (0.62198), and T_v (K) is the virtual dew point temperature.

$$\begin{array}{}\text{(5)}& T_{\mathrm{v}}=T_{\mathrm{dew}}(\mathrm{1}+\left({\displaystyle \frac{\mathrm{1}}{R}}-\mathrm{1.0}\right)q_{\mathrm{0}})\end{array}$$

The height at which the air becomes saturated z is

$$\begin{array}{}\text{(6)}& z={\displaystyle \frac{T_{\mathrm{0}}-T_{\mathrm{dew}}}{{\mathit{\gamma }}_{\mathrm{dry}}}}.\end{array}$$

The elevated temperature following the moist adiabat is then

$$\begin{array}{}\text{(7)}& T_z=T_{\mathrm{dew}}-\left(\mathrm{\Delta }z-z\right){\mathit{\gamma }}_{\mathrm{moist}}.\end{array}$$

Additionally, when T_z<T_dew, 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 the Goff–Gratch formula (Bakan and Hinzpeter, 1987). The adjusted humidity is then used in the surface exchange calculation.

2.1.2 Longwave radiation

Downward longwave radiation is adjusted by assuming the atmosphere behaves as a black body using Stefan–Boltzmann's law. The radiative air temperature at the surface T_rad,0 is calculated using the downward longwave radiation provided by the forcing data LW_↓z0

$$\begin{array}{}\text{(8)}& T_{\mathrm{rad},\mathrm{0}}={\left({\displaystyle \frac{{\mathrm{LW}}_{\downarrow z\mathrm{0}}}{\mathit{\sigma }}}\right)}^{\frac{\mathrm{1}}{\mathrm{4}}},\end{array}$$

where σ is the Stefan–Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴). The radiative temperature at height is then adjusted:

$$\begin{array}{}\text{(9)}& T_{\mathrm{rad},z}=T_{\mathrm{rad},\mathrm{0}}+T_z-T_{\mathrm{0}},\end{array}$$

where T₀ is the grid box mean temperature from the forcing data and T_z is the elevated air temperature. This is used to calculate the downward longwave radiation LW_↓z at height

$$\begin{array}{}\text{(10)}& {\mathrm{LW}}_{\downarrow z}=\mathit{\sigma }{T}_{\mathrm{rad},z}^{\mathrm{4}}.\end{array}$$

An additional correction is made to ensure that the grid box mean downward longwave radiation is preserved:

$$\begin{array}{}\text{(11)}& {\mathrm{LW}}_{\downarrow z}={\mathrm{LW}}_{\downarrow z}-\sum _{z=\mathrm{1}}^n{\mathrm{LW}}_{\downarrow z}\cdot \mathrm{frac}\left(z\right){\displaystyle \frac{\left|{\mathrm{LW}}_{\downarrow z}\right|}{\sum _{i=\mathrm{1}}^z\left|{\mathrm{LW}}_{\downarrow z}\right|\cdot \mathrm{frac}\left(z\right)}},\end{array}$$

where frac is the tile fraction.

2.1.3 Precipitation

To account for orographic precipitation, large-scale and convective rainfall and snowfall are adjusted for elevation using an annual mean precipitation gradient (%/100 m):

$$\begin{array}{}\text{(12)}& P_z=P_{\mathrm{0}}+P_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{precip}}\left(z-z_{\mathrm{0}}\right),\end{array}$$

where P₀ is the surface precipitation, γ_precip is the precipitation gradient, and z₀ is the grid box mean elevation. Rainfall is also converted to snowfall when the elevated air temperature T_z is less than the melting temperature (0 ^∘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 obtain sufficient accumulation in the mass balance compared to observations. The consequence of this is that the grid box mean precipitation is no longer conserved. We tested scaling the precipitation in a way that conserves the grid box mean by reducing the precipitation near the surface and increasing it at height, but this did not yield enough precipitation to obtain a good agreement with the mass balance observations. If the model is being used to simulate river discharge in glaciated catchments, then the precipitation lapse rate could be used as a parameter to calibrate the discharge.

2.1.4 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 effects of katabatic winds on the mass balance when trying to model glacier melt, particularly at lower elevations where the katabatic wind speed is highest.

To explicitly model katabatic winds would require knowledge of the grid box 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

$$\begin{array}{}\text{(13)}& u_z=u_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{wind}},\end{array}$$

where γ_wind is a wind speed scale factor and u₀ 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 obtain 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 longwave radiation, albedo, temperature, wind speed, and roughness length at five heights between 2205 and 3325 m a.s.l. on the glacier. Table S6 in the Supplement 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 1 order of magnitude and the modelled wind speed is 4 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. The surface exchange coefficient, which is used to calculate the sensible heat flux, is a function of the wind speed in the model.

2.2 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). In this scheme the change in albedo as snow ages is related to the growth of the snow grain size, which is a 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 parameterisation which was developed for implementation in the 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 10 cm of the snowpack (ρ_surface) is greater than the firn density (550 kg m⁻³) and the original snow aging scheme is used when (ρ_surface) is less than the firn density.

$$\begin{array}{}\text{(14)}& {\mathit{\alpha }}_{\mathit{\lambda }}={\mathit{\alpha }}_{\mathit{\lambda },\mathrm{ice}}+\left({\mathit{\rho }}_{\mathrm{surface}}-{\mathit{\rho }}_{\mathrm{ice}}\right)\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathit{\lambda },\mathrm{snow}}-{\mathit{\alpha }}_{\mathit{\lambda },\mathrm{ice}}}{{\mathit{\rho }}_{\mathrm{snow}}-{\mathit{\rho }}_{\mathrm{ice}}}}\right)\end{array}$$

α_λ,snow is the maximum albedo of fresh snow, α_λ,ice is the albedo of melting ice, ρ_snow is the density of fresh snow (250 kg m⁻³), and ρ_ice is the density of ice (917 kg m⁻³). The albedo scaling is calculated separately in two radiation bands: visible (VIS) wavelengths λ=0.3–0.7 µm and near-infrared (NIR) wavelengths λ=0.7–5.0 µm. 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 Sect. 3).

2.3 Initialisation

The model requires initial conditions for (1) the snowpack properties and (2) glaciated and unglaciated elevated tile fractions within a grid box. The location of glacier grid points, the initial tile fraction, and the present-day ice mass are 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 twenty-first century. A new feature of the RGI6 is 0.5^∘ gridded glacier volume and area datasets, produced at 50 m 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 and Farinotti, 2012). The area and volume of individual glaciers have been aggregated onto 0.5^∘ grid boxes. We bin the 50 m area and volume into elevations bands varying from 0 to 9000 m in increments of 250 m to match the elevation bands prescribed in JULES.

2.3.1 Initial tile fraction

The elevated glaciated fraction is

$$\begin{array}{}\text{(15)}& {\mathrm{frac}}_{\mathrm{ice}\left(n\right)}={\displaystyle \frac{\mathrm{RGI}\mathit{\_}\mathrm{area}\left(n\right)}{\mathrm{gridbox}\mathit{\_}\mathrm{area}\left(n\right)}},\end{array}$$

where RGI_area is the area (km²) at height from the RGI6, n is the tile elevation, and gridbox_area (km²) is the area of the grid box. In this configuration of the model, any area that is not glaciated is set to a single unglaciated tile fraction (frac_rock) with a grid box mean elevation. It is possible to have an unglaciated tile fraction at every elevation band, but since the glaciated tile fractions do not grow or shrink, we reduce our computation cost by simply putting any unglaciated area into a single tile fraction.

$$\begin{array}{}\text{(16)}& {\mathrm{frac}}_{\mathrm{rock}}=\mathrm{1}-\sum _{n=\mathrm{1}}^{n=n\mathrm{Bands}}{\mathrm{frac}}_{\mathrm{ice}}\left(n\right)\end{array}$$

nBands =37 is the number of elevation bands.

3.1 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 WGMS (2017). Calibrating mass balance against in situ observations is a technique which has been used by other glacier modelling studies (Radić and Hock, 2011; Giesen and Oerlemans, 2013; Marzeion et al., 2012). For the calibration, annual elevation-band mass balance observations are used because there are data available for 16 of the 18 RGI6 regions. For validation, winter and summer elevation-band mass balance is used because there are fewer data available.

The tuneable parameters for mass balance are VIS snow albedo (α_vis,snow), VIS melting ice albedo (α_vis,ice), NIR snow albedo (α_nir,snow), NIR melting ice albedo (α_nir,ice), orographic precipitation gradient (γ_precip), temperature lapse rate (γ_temp), and wind speed scaling factor (γ_wind).

Random parameter combinations are selected using Latin hypercube sampling (McKay et al., 1979) among 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 VIS to NIR albedo using values 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 and 0.7 for which an upper bound value comes from observations of very clean snow with few impurities in the Antarctic (Hudson et al., 2006). The lower bound represents contaminated 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). VIS 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 and 0.5 for which the upper bound comes from spectral albedo observations made in Antarctica (Reijmer et al., 2001). We use a very low minimum albedo for ice in the VIS and NIR wavelengths (0.1), in order to capture the low reflectance of melting ice.

Table 1Tuneable parameters for mass balance calculation and their ranges from the literature.

Download Print Version | Download XLSX

The temperature lapse rate is tuned between values of 4.0 and 10 ^∘C km⁻¹ for which the upper limit is determined from physically realistic bounds and the lower limit is from observations based at glaciers in the 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 of 1–4 to account for an increase in wind speed with height and for 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 Supplement contains the wind speed observations on the Pasterze glacier. The maximum observed wind speed was 4.6 m s⁻¹ (at 2420 m a.s.l.) while the WATCH–ERA Interim dataset (WFDEI) (Weedon et al., 2014) surface wind speed for the same time period was 1.1 m s⁻¹, indicating a scaling factor of approximately 4.

The orographic precipitation gradient γ_precip is tuned between 5 and 25 % per 100 m. This parameter is poorly constrained by observations; therefore a large tuneable range is sampled. Tawde et al. (2016) estimated a precipitation gradient of 19 % per 100 m for 12 glaciers in the western Himalayas using a combination of remote sensing and in situ meteorological observations of precipitation. Observations show that the precipitation gradient can be as high as 25 % per 100 m for glaciers in Svalbard (Bruland and Hagen, 2002) while glacier–hydrological modelling studies have used much smaller values of 4.3 % per 100 m (Sorg et al., 2014) and 3 % per 100 m (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 humidity, rainfall, snowfall, and wind speed from the WFDEI dataset (Weedon et al., 2014). To reduce the computation time, only grid points at which glacier ice is present are modelled. An ensemble of 198 calibration experiments are run. For each simulation the model is spun up for 10 years and the elevation-dependent mass balance is compared to observations at 149 field sites over the years 1979–2014.

The elevation-dependent mass balance observations come from stake measurements taken every year at different heights along the glaciers. Many of the mass balance observations in the World Glacier Monitoring Service (WGMS, 2017) are supplied without observational dates. In this case, we assume the mass balance year starts on 1 October and ends on 30 September with the summer commencing on 1 May. Dates in the Southern Hemisphere are shifted by 6 months. The observations are grouped according to standardised 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.

Figure 2The location of mass balance profile observation glaciers from the World Glacier Monitoring Service and the Randolph Glacier Inventory regions (version 6.0).

Download

Figure 3Modelled annual elevation-dependent specific mass balance against observations from the WGMS. The modelled mass balance is simulated on a 0.5^∘ grid resolution at 250 m elevation bands and the observations are for individual glaciers at elevation levels specific to each glacier. The observed mass balance is interpolated onto the JULES elevation bands. If only a single observation exists, then mass balance for the nearest JULES elevation band is used. The number of glaciers is shown in the top left-hand corner and the number of observation points in brackets. In central Europe mass balance for the Maladeta glacier in the Pyrenees is shown in black circles.

Download

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 16 regions have a negative bias in the annual mass balance. Notably Svalbard, southern Andes, and New Zealand underestimate mass balance by 1 m w.e. yr⁻¹. The negative bias is also seen in the summer and winter mass balance and discussed in Sect. 3.2. The model performs particularly poorly for the low-latitude region, which has a large RMSE (3.02 m w.e. yr⁻¹). This region 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 Climatic Research Unit (CRU) data (Mitchell and Jones, 2005). 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 includes sublimation, so it is possible the WFDEI data over tropical glaciers are too warm. The WFDEI data are based on the ERA-Interim reanalysis in which air temperature has been constrained using CRU data. The CRU data comprise 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. In central Europe some of the poor correlations with observations are caused by the Maladeta glacier in the Pyrenees (Fig. 3), which is a small glacier with an area of 0.52 km² (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 m of water equivalent per year.

Table 2Best parameter sets for each RGI6 region. The regions are ranked from the lowest to the highest RMSE. There are no observed profiles for Iceland and the Russian Arctic, so the global mean parameter values are used (bold) for the future simulations.

Download Print Version | Download XLSX

Table 3Root-mean-square error (RMSE), correlation coefficient (r), Nash–Sutcliffe efficiency coefficient (NS), mean bias (Bias), and the number of elevation-band mass balance observations (No. of obs) for RGI6 regions. The regions are ranked from the lowest to the highest RMSE.

Download Print Version | Download XLSX

3.2 Model validation

The calibrated mass balance is validated against summer and winter elevation-band specific mass balance for each region where data are available (Fig. 4). For all regions, except Scandinavia in the summer, negative Nash–Sutcliffe 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. This means that future projections of volume loss presented in Sect. 4.2 might be overestimated.

Figure 4Comparison 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.

Download

Table 4Winter (bold) and summer number of elevation-band mass balance observations (No. of obs), root-mean-square error (RMSE), correlation coefficient (r), Nash–Sutcliffe efficiency coefficient (NS), and mean bias (Bias).

Download Print Version | Download XLSX

The reason for the negative bias is because the model underestimates the precipitation and therefore the accumulation part of the mass balance is underestimated. This is because our approach to correcting the coarse-scale gridded precipitation for orographic effects is simple. We use a single precipitation gradient for each RGI6 region and do not apply a bias correction. A bias correction is often recommended because precipitation is underestimated in coarse-resolution datasets. Gauging observations are sparse in high-mountain regions and snowfall observations can be susceptible to undercatch by 20 %–50 % (Rasmussen et al., 2012). Our precipitation rates are generally too low because we do not bias-correct the precipitation.

Other studies use a bias correction that varies regionally (Radić and Hock, 2011; Radić et al., 2014; Bliss et al., 2014). In those studies, the precipitation at the top of the glacier was estimated using a bias correction factor k_p. The decrease in precipitation from the top of the glacier to the snout was calculated using a precipitation gradient. To account for the fact that the mass balance of maritime and continental glaciers responds differently to precipitation changes, k_p was related to a continentality index. Our motivation for using a single precipitation gradient for each RGI6 region and no bias correction was to test the simplest approach first; however the resulting biases suggest that this approach could be improved.

The impact of underestimating the precipitation is that we simulate negative mass balance in winter at some observational sites (Figs. 5a and 4). To demonstrate this, we compare the mass balance components for two glaciers: the Leviy Aktru in the Russian Altai Mountains, which has negative mass balance in the winter, and Kozelsky glacier in northeastern Russia, which has no negative mass balance in the winter (See Fig. S9). Both glaciers are in the north Asia RGI6 region, so they have the same tuned parameters for mass balance. The simulated winter accumulation rates are much lower at Leviy Aktru glacier than Kozelsky glacier, leading to negative mass balance at the lowest three model levels below 2750 m.

Figure 5Simulated and observed elevation-dependent winter mass balance when grid boxes with a glacier area of less than 100, 300, and 500 km² are excluded. The colour identifies the RGI6 regions shown in Fig. 2. The RMSE, correlation coefficient, and number of glaciers are listed.

Download

Figure 6Simulated and observed elevation-dependent summer mass balance when grid boxes with a glacier area of less than 100, 300, and 500 km² are excluded. The colour identifies the RGI6 regions shown in Fig. 2. The RMSE, correlation coefficient, and number of glaciers are listed.

Download

The simplistic treatment of the precipitation lapse rate also leads to instances in which the model simulates positive mass balance in the summer at some locations (Figs. 6a and 4). We show the summer mass balance components for the same two glaciers in Fig. S10. Positive mass balance is simulated at Kozelsky glacier because accumulation exceeds the melting. This suggests that the precipitation gradient (19 % per 100 m for north Asia) is overly steep in the summer at this location.

Another reason we underestimate the accumulation is due to the partitioning of rain and snow based on an air temperature threshold of 0 ^∘C. The 0 ^∘C 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 (W_k) and ice content (I_k) in each layer (k)

$$\begin{array}{}\text{(17)}& C_k=I_k{C}_{\mathrm{ice}}+W_k{C}_{\mathrm{water}},\end{array}$$

where C_ice=2100 JK⁻¹ kg⁻¹ and C_water=4100 JK⁻¹ kg⁻¹. The liquid water content is limited by the available pore space in the 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.5 ^∘C (Huss and Hock, 2015; Giesen and Oerlemans, 2012), 2 ^∘C (Hirabayashi et al., 2010), and 3 ^∘C (Marzeion et al., 2012). An improved approach would use the wet-bulb temperature to partition rain and snow, which would include the effects of humidity on temperature. Alternatively, a spatially varying threshold based on precipitation observations could be used. Jennings et al. (2018) showed, by analysing precipitation observations, that the temperature threshold varies spatially and is generally higher for continental climates than maritime climates.

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, performs better in the winter when snow is accumulating than in summer when there is melting. Figures 5b–d and 6b–d show the winter and summer mass balances for all observation sites when area thresholds of 100, 300, and 500 km² are applied to the validation. There is an improvement in the simulated summer mass balance when the glaciated area increases. This is seen by the improved correlation in Fig. 6d in which the validation is repeated but only grid boxes with a glaciated area greater than 500 km² are considered. This indicates the model is better at simulating summer melting over regions with a large ice extent than over regions with a small glaciated area.

Table 5List of high-end climate change CMIP5 models that are downscaled using HadGEM3-A. The years when the CMIP5 models pass +1.5, +2, and +4 ^∘C global average warming relative to the pre-industrial period are shown.

* No data are available for 2113 because the bias-corrected data end at 2097.

Download Print Version | Download XLSX

4.1 Downscaled climate change projections

Glacier volume projections are made for all regions, excluding Antarctica, for a range of high-end climate change scenarios. This is defined as climate change that exceeds 2 and 4 ^∘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). 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 ^∘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 consider for glaciers in the central and eastern Himalayas, 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 (ISIMIP) (Hempel et al., 2013). This technique uses WATCH forcing data (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 bilinearly interpolated from its native resolution of N216 (∼60 km), onto a 0.5^∘ grid, to match the resolution of the WATCH forcing data, which were used for the bias correction. The daily bias-corrected surface fields from the HadGEM3-A are used to run JULES offline to calculate future glacier volume changes. The bias correction was only applied to data up until the year 2097, which means the glacier projections terminate at this year. A flow chart of the experimental setup is shown in Fig. 7. The HadGEM3-A climate data were generated and bias-corrected for the High-End cLimate Impact and eXtremes (HELIX) project.

Figure 7Flow chart showing the experimental setup to calculate future glacier volume. * The bias correction method is described by Hempel et al. (2013).

Download

Table 6Percentage ice volume loss, relative to the initial volume (ΔV), and ice loss in millimetres of sea level equivalent (SLE) for the end of the century (2097). Percentage volume losses are shown for low, medium, and high elevation ranges as well as for all elevations. The data show the multi-model mean ± 1 standard deviation. The conversion of volume to SLE assumes an ocean area of 3.618×10⁸ km². The initial area and volume from the Randolph Glacier Inventory version 6 are listed in columns 1 and 2.

Download Print Version | Download XLSX

4.2 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 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 the end of the run (2097) to the initial volume calculated from the RGI6. Regional volume changes expressed in percent for low (0–2000 m), medium (2250–4000 m), high (4250–9000 m), and all elevation ranges (0–9000 m) are listed in Table 6. The total volume loss over all elevation ranges is also listed in millimetres of sea level equivalent in Table 6 and plotted in Fig. 10. Maps of the percentage volume change at the end of the century, relative to the initial volume, are contained in the Supplement in Figs. S1–S7.

Figure 8Regional glacier volume projections using the HadGEM3-A ensemble of high-end climate change scenarios.

Download

A substantial reduction in glacier volume is projected for all regions (Fig. 8). Global glacier volume is projected to decrease by 64±5 % by end of the century, for which the value corresponds to the multi-model mean ± 1 standard deviation. The regions which lose more than 75 % of their volume by the end of the century are Alaska ($-\mathrm{89}\pm \mathrm{2}$ %), western Canada and the US ($-\mathrm{100}\pm \mathrm{0}$ %), Iceland ($-\mathrm{98}\pm \mathrm{3}$ %), Scandinavia ($-\mathrm{98}\pm \mathrm{3}$ %), the Russian Arctic ($-\mathrm{79}\pm \mathrm{10}$ %), central Europe ($-\mathrm{99}\pm \mathrm{0}$ %), Caucasus ($-\mathrm{100}\pm \mathrm{0}$ %), central Asia ($-\mathrm{80}\pm \mathrm{7}$ %), South Asia west ($-\mathrm{98}\pm \mathrm{1}$ %), South Asia east ($-\mathrm{95}\pm \mathrm{2}$ %), low latitudes (100±0 %), southern Andes ($-\mathrm{98}\pm \mathrm{1}$ %), and New Zealand ($-\mathrm{88}\pm \mathrm{5}$ %). The HadGEM3-A forcing data show these regions experience a strong warming. In most regions this is combined with a reduction in snowfall relative to the present day, which drives the mass loss (Fig. 9). Regions most resilient to volume losses are Greenland ($-\mathrm{31}\pm \mathrm{5}$ %) and Arctic Canada north ($-\mathrm{47}\pm \mathrm{3}$ %). 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 glacier models (Marzeion et al., 2012; Huss and Hock, 2015). Our model overestimates the melting of these glaciers for the calibration period (Fig. 3), so this result should be treated with a degree of caution. Some of the high-latitude regions, particularly Alaska, western Canada and the US, Svalbard, and north Asia, experience very large volume increases at their upper elevation ranges. This would be reduced if the model included glacier dynamics, because ice would be transported from higher elevations to lower elevations. The ensemble mean global sea level equivalent contribution is 215.2±21.3 mm. The largest contributors to sea level rise are Alaska (44.6±1.1 mm), Arctic Canada north and south (34.9±3.0 mm), the Russian Arctic (33.3±4.8 mm), Greenland (20.1±4.4), high-mountain Asia (combined central Asia and South Asia east and west), (18.0±0.8 mm), southern Andes (14.4±0.1 mm), and Svalbard (17.0±4.6 mm). These are the regions which have been observed by the Gravity Recovery and Climate Experiment (GRACE) satellite to have lost the most mass in the recent years (Gardner et al., 2013).

Figure 9Regional temperature and snowfall changes relative to the present day (2011–2015) from the HadGEM3-A ensemble over glaciated grid points. The ensemble mean is shown in the solid line and the range of model projections are shown in the shaded regions.

Download

To investigate which parts of the energy balance are driving the future melt rates, we show the energy balance components averaged over all regions and all elevation levels in Fig. 11. Future melting is caused by a positive net radiation of approximately 30 W m⁻² that is sustained throughout the century. This is comprised of 18 W m⁻² net shortwave radiation, 3 W m⁻² net longwave radiation, 5 W m⁻² latent heat flux, and 4 W m⁻² sensible heat flux. The largest component of the radiation for melting comes from the net shortwave radiation. The upward shortwave radiation comprises direct and diffuse components in the VIS and NIR wavelengths. The VIS albedo decreases because melting causes the ice surface to darken. In contrast, the NIR albedo increases because the ice is heating up, emitting radiation in the infrared part of the spectrum.

Figure 10(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.

Download

Figure 11Ensemble mean energy balance components averaged over all glaciated regions and all elevation bands when the model is forced with HadGEM3-A data.

Download

The downward and upward longwave radiation are increasing in the future; however, the net longwave radiation contribution to the melting is small. The downward longwave radiation increases because of the T⁴ relationship with air temperature, whereas the upward longwave radiation increases because the glacier surface is warming. The latent heat flux from refreezing of meltwater and the sensible heat from surface warming are also small components of the net radiation balance.

4.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 elevation ranges for the period 1980–2000. Sublimation is excluded because its contribution to mass balance is relatively 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 meltwater 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 increase 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.e. yr⁻¹) and western Canada and the US (7.3 m w.e. yr⁻¹) between 4250 and 9000 m and the largest melt rates are found in the Caucasus and Middle East (−7.4 m w.e. yr⁻¹) and the low latitudes (−7.6 m w.e. yr⁻¹).

Figure 12Modelled 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 a positive sign and sublimation is excluded because its contribution is very small. Mass balance components are averaged over low (0–2000 m), medium (2250–4000 m), and high (4250–9000 m) elevation ranges.

Download

Figure 13Global mass balance components for three elevation ranges. The historical period is calculated using the WFDEI 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.

Download

Figure 13 shows how the global annual mass balance components vary with time for low, medium, and high elevation ranges.

There is a reduction in accumulation and refreezing at all elevation ranges towards the end of the century. Melt rates decrease at medium and high elevation ranges because glacier mass is lost at these altitudes; therefore less ice is available to melt (see Fig. S8 for the future cumulated mass balances as a function of height). Melt rates are constant at the low elevation ranges because there remain substantial quantities of ice available to melt at the end of the century in Greenland, Arctic Canada north and south, Svalbard, and the Russian Arctic. At high elevations mass balance is reduced from −2.2 m w.e. yr⁻¹ (−177 Gt yr⁻¹) during the historical period (1980–2000) to −0.35 m w.e. yr⁻¹ (−28 Gt yr⁻¹) by the end of the century (2080–2097). Similarly, for the medium elevation ranges mass balance decreases from −0.56 m w.e. yr⁻¹ (−26 Gt yr⁻¹) to −0.24 m w.e. yr⁻¹ (−11 Gt yr⁻¹).

4.4 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 sea surface temperature and sea ice concentration boundary conditions, provided by the CMIP5 models, on the HadGEM3-A 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, 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 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 change-point detection algorithm called findchangepts (Rebecca et al., 2012; Lavielle, 2005) from the MATLAB signal processing toolbox. The algorithm is run to find where the mean of the top 10 experiments (excluding the optimal experiment) changes the most significantly. For each RGI6 region the step changes in the RMSE are shown in Fig. 14.

Figure 14Calibration experiments ranked 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.

Download

JULES is rerun for each of the downscaled CMIP5 experiments and for each parameter set that is defined as performing equally well (see Table S1 in the Supplement for a list of the parameter sets that perform equally well). The volume losses expressed in millimetres of sea level equivalent are shown in Fig. 15. The effects of the parametric uncertainty on the volume losses vary regionally, with the largest impact found for central Europe and Greenland. Regional volume losses 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 of sea level equivalent volume loss by the end of the century.

Figure 15Regional 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.

Download

Table 7Regional ensemble mean, minimum, and maximum volume losses for 2097 in sea level equivalent (mm) when the present-day 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.

Download Print Version | Download XLSX

Another way to explore the uncertainty in the volume projections caused by the calibration procedure is to use different performance metrics to identify the best parameter sets. In addition to using RMSE, we calculate the best parameter sets by (1) minimising the absolute value of the bias and (2) maximising 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 Supplement). For 12 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 above, this could be due to the fact the grid box mean WFDEI precipitation was not bias-corrected. Glacier volume projections are generated by repeating the simulations using these two additional performance metrics to identify the 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 of sea level equivalent by the end of the century.

Figure 16Multi-model mean (black line) and ensemble spread (shaded) global volume loss in sea level equivalent. Panel (a) shows the volume loss when optimum parameter sets are selected by minimising the RMSE and (b) is volume loss when optimum parameter sets are selected using additional performance metrics (minimising RMSE, minimising the bias, and maximising the correlation coefficient).

Download

4.5 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 (Huss and Hock, 2015; Radić 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) listed regional percentage volume change and sea level equivalent values in their study while Radić et al. (2014) listed sea level equivalent values only (see the comparison of Tables S4 and S5 in the Supplement).

Our end-of-the-century percentage volume losses compare reasonably well to Huss and Hock (2015) for central Europe, Caucasus, South Asia east, Scandinavia, the Russian Arctic, western Canada and the 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 south. The uncertainty in our percentage volume losses is smaller than that of 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±20 mm, which is higher than 188 mm (Radić et al., 2014) and 136±23 mm (Huss and Hock, 2015) caused mainly by greater contributions from Alaska, southern Andes, and the Russian Arctic. These three regions are discussed in turn.

For the southern Andes our estimates are approximately double (14.4 mm) those of the other studies (5.8 mm, Huss and Hock, 2015; 8.5 mm, Radić et al., 2014). This region has the largest negative bias in the calibrated present-day mass balance (−2.87 m w.e. yr⁻¹; 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 grid box is assumed to have the same regional mass balance bias. The bias-corrected volume losses are listed in Table S5 in the Supplement. For the southern Andes, the volume losses are much closer to the other studies (7.6 mm) 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 obtain a tuned parameter set (global mean parameters were used instead). In Alaska the bias in annual mass balance is small (0.06 m w.e. yr⁻¹), 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±20 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 sea level equivalent contributions are higher than the other studies because the climate forcing data are different. The HadGEM3-A model uses boundary conditions from a subset of RCP8.5 CMIP5 models with the highest warming levels. Furthermore, the HadGEM3-A data have a higher resolution (approximately 60 km) than the CMIP5 data, which were 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.