2.1. Morteratsch glacier complex

The 16 km² glacier complex is located in southeastern Switzerland (tongue at 46.43°N, 9.93° E in 2012) and consists of two distinct glaciers: the main Vadret da Morteratsch (with a length of 6 km) and its tributary, Vadret Pers (~4.5 km) (Fig. 1). At present, the glacier front is at an elevation of ~2100 m, while the highest parts are ~2000 m higher and correspond to mountain tops such as Piz Bernina (4049 m), Piz Zupo (3996 m) and Piz Palü (3905 m). These peaks prevent direct insolation on the glacier by shadowing for a large part of the year and therefore affect the surface mass balance.

Fig. 1. Annual surface mass balance averaged over the period 2001–10. The map uses the Swiss CH1903 coordinate system. The thick black line delineates the Morteratsch glacier complex; the thin black lines are surface-elevation contours at 200 m intervals (indicated from 2200 to 3800 m). Note that only ice that flows into the glacier complex is considered in the mask, and that isolated small glaciers, such as Vadret da la Fortezza (situated between Vadret Pers (P) and Vadret da Morteratsch (M)), are not taken into account. The highest surrounding peaks are Piz Bernina (PB, 4049 m), Piz Zupo (PZ, 3996 m) and Piz Palü (PP, 3905 m). The inset in the upper right corner shows the location of the glacier in Switzerland.

2.2. Field data

Starting in 2001, we have collected an elaborate dataset on field observations. The ice thickness was measured with different ground-penetrating radar (GPR) systems and subsequently extrapolated over the entire glacier assuming plastic flow along central flowlines. Based on this, we estimate that the glacier volume in 2001 was 1.14 ± 0.28 km³ (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). Since 2001, we have measured the annual mass balance with ablation stakes, resulting in a dataset of 186 individual measurements so far, which have also been used to deduce the surface velocity of the glacier (error within the order of 1 m a⁻¹). We also performed accumulation measurements with aluminium poles. The absolute movement of all stakes is known from GPS measurements (error within the order of 1 m in the three dimensions). For sites near the ice front where horizontal displacements are virtually nil, this directly provides the thinning rate.

2.3. Ice-dynamic model

A 3-D isothermal ice-dynamic model (Reference Fürst, Rybak, Goelzer, De Smedt, De Groen and HuybrechtsFürst and others, 2011) with Blatter/Pattyn-type equations (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003) is used, for which the rate factor for temperate ice A ₀ and sliding parameter A _sl were tuned to match observed surface velocities (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). Nye’s generalization of Glen’s flow law is used to describe the internal deformation (Reference GlenGlen, 1955; Reference NyeNye, 1957):

where τ_ij are the deviatoric stresses, n is the power-law exponent (equal to 3 in this study) and is a small offset of 10^–30 to prevent infinite effective viscosity (η) when the strain rate is zero (following Reference PattynPattyn, 2003). The effective stress is determined by the second invariant of the strain-rate tensor:

for which the components of the strain rate tensor ε _ij are defined as

and u_i are the 3-D components of the velocity vector.

Equations (1–3) are simplified to obtain an LMLa (multi-layer longitudinal stresses) higher-order model (Reference HindmarshHindmarsh, 2004; Reference Fürst, Rybak, Goelzer, De Smedt, De Groen and HuybrechtsFürst and others, 2011). Assuming local cryostatic equilibrium, bridging effects are neglected in the force balance, and horizontal gradients in the vertical velocity field are assumed to be negligible in comparison with vertical gradients of horizontal velocity. Under these approximations, the horizontal and vertical velocity fields are decoupled, and the latter can be determined via mass conservation. The basal drag τ _b corresponds to the sum of all basal resistive forces and is calculated following the higher-order approximation (Reference WeertmanWeertman, 1964; Reference Van der Veen and WhillansVan der Veen and Whillans, 1989):

for which b is the bedrock elevation and τ_ij are the components of the stress tensor at the base. The best fit between observed and modelled velocities (root-mean-square error (RMSE) of 15.0 m a⁻¹) was obtained for a rate factor A ₀ of 1.6 × 10⁻¹⁶ Pa⁻³ a⁻¹ and a sliding parameter A _sl of 12 × 10⁻¹⁶ m⁷ N⁻³ a⁻¹ (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). The highest flow velocities are reached for Vadret da Morteratsch and go up to 125 m a⁻¹ (see Fig. 7 in Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). For Vadret Pers, both the rate factor A ₀ and sliding parameter A _sl were doubled for a better match between observed and modelled velocities and to avoid Vadret Pers becoming too thick in transient simulations. Flow-law parameters are kept constant in time. The model is run on a 25 m horizontal resolution, and in the vertical the grid uses 21 non-equally spaced layers. The flow model is coupled to a surface energy-balance model that determines the spatio-temporal distribution of accumulation and ablation.

The continuity equation is the link between ice dynamics and the surface mass balance:

where H is the ice thickness, ū is the vertically averaged horizontal velocity vector, ▽ is the 2-D divergence operator and M is the surface mass balance. The higher-order velocity field and changes in ice thickness are calculated with a weekly time step (Δt = 0.02 years), while the mass balance is updated annually. More details on parameter choices, boundary conditions and numerical implementation are provided by Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others (2013) and Reference Fürst, Rybak, Goelzer, De Smedt, De Groen and HuybrechtsFürst and others (2011).

2.4. Mass-balance model and climatic input

To compute the surface mass balance we use a 2-D energy-balance model (Reference Oerlemans, AA and Lisse OerlemansOerlemans, 2001) that has been calibrated for the Morteratsch glacier complex (Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others, 2009). The model is based on an equation for the cumulative mass balance B (m w.e. a⁻¹) and accounts for the shading effect of the surrounding mountains, the daily temperature cycle, the timing of precipitation events, and seasonal changes in snow thickness and albedo:

Here L _m is the latent heat of melting (3.34 × 10⁵ J kg⁻¹), ρ _w is the water density (1000 kg m⁻³) and P _solid(t) is the solid precipitation. The time step Δt is 1 hour.

The net energy flux ψ (W m⁻²) at the surface is represented by

where the first term on the right-hand side is the net shortwave radiation received at the glacier surface, with τ the transmissivity of the atmosphere and α the surface albedo. Q (W m⁻²) is the incoming solar radiation, which consists of a direct and an indirect (diffuse) term and is corrected for the slope and orientation of the glacier surface (Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others, 2009). The term c ₀ + c ₁ . T _air represents the sum of the net longwave radiation and the turbulent heat exchange (Reference Oerlemans, AA and Lisse OerlemansOerlemans, 2001). Here, c ₀ = –45 W m⁻², c ₁ = 12 W m⁻² °C⁻¹, τ = 0.45 (Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others, 2009) and T _air is the air temperature (°C). The parameter values, in particular for τ, were obtained from a best fit of the energy-balance model with our observed mass-balance measurements, and are further informed by detailed observations from an automatic weather station placed in the lower ablation area (personal communication from J. Oerlemans, 2008). The fraction of solid precipitation depends on the local air temperature. Total precipitation is upscaled from nearby weather stations using a fixed altitudinal gradient, and accounts for the timing and intensity of precipitation events informed by a daily precipitation series from a nearby valley station for the period 1981–2006. As discussed in Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others (2009), available data do not allow a further differentiation of snow accumulation with respect to wind exposure, for example, making the precipitation treatment probably the largest source of uncertainty in the calculation of the net mass balance. The cumulative mass balance B after a full balance year (1 October to 30 September) gives the annual specific mass balance, M, needed in Eqn (6).

Meteorological data from surrounding weather stations were used to calibrate the model to ablation and accumulation measurements on both Morteratsch and Pers glaciers for the period 2001–06 (Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others, 2009). Most measurements are available for the ablation zone, while for the accumulation zone we rely on a handful of measurements from aluminium poles, snow pits and a shallow ice core that was drilled by a Swiss group under Piz Zupo (Reference Sodemann, Palmer, Schwierz, Schwikowski and WernliSodemann and others, 2006). The climatic input to drive the mass-balance model consists of a monthly temperature record from Segl Maria (46.44 ° N, 9.77 ° E, 1798 m a.s.l.) since 1864 and monthly precipitation measurements from Samedan (46.53 ° N, 9.88 ° E, 1705 m a.s.l.) extending back to 1861 (Fig. 2). These data originate from the HISTALP database (Reference AuerAuer and others, 2007) and were homogenized to reduce errors such as the early instrumental warm-bias (Reference Frank, Büntgen, Böhm, Maugeri and EsperFrank and others, 2007; Reference Böhm, Jones, Hiebl, Frank, Brunetti and MaugeriBöhm and others, 2010). The modelled annual mass balance field-averaged over the decade 2001–10, using the topography from a digital elevation model (DEM) of 2001 (DHM25 from SwissTopo; uncertainty 3–6 m), is shown in Figure 1. For this period, the Morteratsch glacier complex has a mean specific mass balance of –0.93 m w.e. a⁻¹ and an accumulation-area ratio (AAR) of 44%, clearly indicative of a glacier out of equilibrium with the imposed climate due to a disproportionally large snout characterized by high ablation. The associated average equilibrium-line altitude (ELA) over this period is ~3020 m for Vadret da Morteratsch, while for Vadret Pers it is ~3080 m. The higher ELA and generally more negative surface mass balance for Vadret Pers are due to glacier orientation. For a large part of the day, most of Vadret Pers is directly exposed to the sun, whereas Vadret da Morteratsch is oriented to the north and receives more shading from the high surrounding peaks (especially from Piz Zupo before noon and from Piz Bernina in the afternoon).

Fig. 2. Annual precipitation measured in Samedan (1861–2010) and average summer temperature (May–September) as recorded in Segl Maria (1864–2010), corrected for early instrumental warm-bias. The homogenized data come from the HISTALP database. The blue lines represent annual values, and the red line represents the 20 year running mean.

2.5. Specific model set-up

Several rock outcrops exist around the ELA and within the accumulation area of the glacier. Most of them are persistent features since the LIA, as identified on topographical maps (1: 50 000) from 1851 (Reference Maisch, Burga and FitzeMaisch and others, 1999) and 1935 (DEM provided by F. Paul). These ice-free areas are very steep, and consequently icefall and wind drift prevent permanent snow or ice cover. Our glacier model does not account for the details of this effect, so these parts are kept ice-free in all simulations. However, it is reasonable to assume that snow falling on these ice-free parts will ultimately end up on the glacier and thus contribute to its mass balance. In the model this is achieved by redistributing the positive mass balance of these regions over an elevation band down to 50 m below the specific ice-free part. An additional issue related to keeping these zones ice-free is that mass is lost when ice flows into these regions, which is the case when the downward surface slope is directed towards the outcrops. In reality this process coincides with ice avalanches that fall on lower parts of the glacier. To take this into account, a similar method to that for the snowfall is used and the mass that flows into the ice-free zones is redistributed over the same elevation bands as for the snow. As the glacier geometry hardly changes in the accumulation zone and to avoid ice flowing over ridges away from the glacier (e.g. to the south, towards Italy), the glacier is not allowed to expand laterally above 3000 m of elevation. These modifications ensure that mass is conserved over the glacier complex at all times.

2.6. Past glacier evolution (1850s to 2010)

An array of observational data is available to reconstruct the retreat history of Morteratsch glacier since the LIA. These data are required for the dynamical calibration of the ice-flow model before any future projections.

3.1. LIA extent of the glacier

The initial condition of the simulation is a glacier in steady state at the LIA extent. Although there is no reason to a priori suspect a steady state around 1864, it is imposed to avoid an unwanted model drift at the start of simulations. Figure 3 shows the corresponding geometry resulting from the flow model. To obtain this glacier extent, the surface mass balance was calculated with average climate conditions for the period 1864–93, which is the first 30 year period for which homogenized HISTALP data from nearby stations are available. However, to obtain the desired glacier length it was necessary to apply an additional mass-balance forcing of +0.5 m w.e. a⁻¹, corresponding to an ELA lowering of ~60 m. This seems reasonable, as the mass balance over the glacier was very likely more positive at the end of the LIA.

Fig. 3. Modelled steady-state thickness field for a glacier with the same length as indicated by the LIA moraine. The contours are spaced 50 m apart. The thick black line indicates the ELA. White areas are either ice-free or contain ice that does not contribute to the flow of the glacier complex.

For the modelled LIA geometry the glacier volume is 2.25 km³ and the thickest ice for central Morteratsch glacier is found at 410 m. The ELA is ~200 m lower than the 2001– 10 average at ~2820 m for Vadret da Morteratsch and 2880 m for Vadret Pers. As a result of the steady-state configuration the mean specific mass balance of the glacier complex, which is the annual surface mass balance averaged over the whole glacier, is equal to zero.

3.2. Evolution of the glacier from 1864 to 2010

From the balance year 1864/65 onwards, the surface mass balance B is calculated from the 2-D energy-balance model driven by the Segl-Maria and Samedan temperature and precipitation records. In the dynamic calibration process, a mass-balance correction δB is iteratively applied to best match the observed evolution of the frontal position as follows:

Here B is a function of location (x, y) and time t, B _SMB is calculated from the energy-balance model and δB (t) is considered uniform over the glacier. The mass-balance correction implicitly accounts for deviations from the initial steady-state assumption, uncertainties in the mass-balance model and errors in the climatic-forcing records. This correction varies in time and is necessary to reproduce the observed retreat. This approach to generate the current glacier imbalance is similar to Reference OerlemansOerlemans (1997) in a study on Nigardsbreen. Norway, and is applicable when the glacier-length record exceeds the characteristic response time (Reference Zuo and OerlemansZuo and Oerlemans, 1997; Reference De Smedt and PattynDe Smedt and Pattyn, 2003).

Figure 4 shows the required mass-balance correction δB (t), which was determined manually in order to reproduce the main frontal retreat trends. It can conveniently be expressed as four piecewise linear sections representing different periods between 1864 and 1960. The initial δB is kept constant at +0.5 m w.e. a⁻¹ until 1905. Without keeping this positive mass-balance correction, the modelled retreat would have been too strong. After 1905, δ B had to decrease to 0.1 m w.e. a⁻¹ by 1930 to capture the accelerated retreat starting in the 1930s. A further positive mass-balance correction of maximally +0.5 m w.e. a⁻¹ centred around 1945 was required to correctly simulate the smaller glacier retreat rates during the 1970s and 1980s. In this paper we can only speculate about the reasons why a surface-mass-balance model driven by HISTALP data alone is not able to produce a correct glacier retreat. Apparently, the glacier’s mass balance is affected by additional features not captured well in the ambient climate records, as noted earlier, for example in a study on Glacier d’Argentière, France (Reference Huybrechts, De Nooze, Decleir and OerlemansHuybrechts and others, 1989), and a recent study by Reference Lüthi, Bauder and FunkLüthi (2014). Possibly this is related to variations in aerosol content (Reference Huss, Funk and OhmuraHuss and others, 2009; Reference Painter, Flanner, Kaser, Marzeion, Van Curen and AbdalatiPainter and others, 2013). The absence of a mass-balance correction after 1960 to correctly simulate the glacier retreat of the last 50 years does, however, lend credibility to the dynamic calibration as performed here. It clearly suggests that the effect of the initial steady-state assumption has faded by 2010 and that the surface-mass-balance model driven by unbiased HISTALP data over recent decades produces a basically correct forcing for the glacier flow model.

Fig. 4. Modelled and observed glacier length evolution of Vadret da Morteratsch between 1878 and 2010. An additional mass-balance correction (shown in blue) is required for a better match of the modelled frontal retreat with the observations and is set to zero after 1960.

Other sources of information have been used to validate the glacier model. Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) reconstructed the evolution of the Morteratsch glacier volume between 1900 and 2008 based on mass-balance modelling and absolute volume differences derived from four DEMs for 1935, 1955, 1985 and 2008. Overall, there is a good agreement between the results of Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) and our volume change reconstructions (Fig. 5), especially after 1960 when the mass-balance correction is no longer applied. In both cases, the same general trends are captured, including two phases of higher volume-loss rates. Interestingly, these two phases (1925–1960 and from 1980 onwards) precede the associated periods of accelerated frontal retreat by ~5–10 years. This is because glacier volume is known to adapt faster to climatic changes than glacier length (e.g. Reference Oerlemans, AA and Lisse OerlemansOerlemans, 2001; Reference Adhikari, Marshall and HuybrechtsAdhikari and others, 2011). This can also be inferred from the evolution of glacier area shown in Figure 5, which mainly reflects length changes at the snout. Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) derived a volume loss of 0.61 km³ between 1935 and 2008, compared with 0.58 km³ from our ice-dynamic model. For the more recent period between 1985 and 2008, the volume losses are respectively 0.305 km³ and 0.316 km³. Reference Fischer, Huss and HoelzleFischer and others (2014), also based on DEM differencing, obtain a volume loss of 0.30 km³ between 1991/92 and 2008/09, while we obtain 0.27 km³ for this period.

Fig. 5. Reconstructed cumulative volume change (relative to 2008) and area change (2 year running mean) for the Morteratsch glacier complex between 1864 and 2010. The reconstruction by Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) is based on mass-balance modelling, constrained by four fixed points derived from DEM differencing.

The recent increase in mass loss and the subsequent accelerated frontal retreat has its counterpart in the long-term mean specific mass balance of the Morteratsch glacier complex. Centred around 1980, the 20 year average annual surface mass balance of the glacier is around –0.21 m w.e. a⁻¹, while only 20 years later, in 2000, this value lowers to –0.82 m w.e. a⁻¹. From the meteorological records shown in Figure 2, it is clear that the recent glacier retreat is primarily driven by higher summer temperatures.

Additionally, the modelled 2001 glacier state is compared with the observed glacier geometry (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). As shown in Figure 6, there is overall quite good agreement between the modelled surface topography and that from the DEM. As the bedrock is kept constant in the modelling, good agreement also applies to the ice thickness. The highest ice thickness of ~350 m is found for central Vadret da Morteratsch in both fields. The maximal ice thickness for Vadret Pers is equally 250 m in both cases. The total ice volume obtained from the transient model run in 2001 (1.41 km²) is, however, situated at the high end of the range inferred from available radar observations (1.14 ± 0.28 km³) (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). This may reflect shortcomings in the assumptions made to generalize the radio-echo sounding lines over the whole glacier to reconstruct ice thickness, as well as in the flow parameters used in the ice-dynamic model. The largest discrepancies between observed and modelled ice thickness occur for the tongue, the confluence area and the accumulation area of Morteratsch glacier. However, the last zone lacks meaningful ice thickness and velocity measurements. One explanation is that the flow model assumption for the ice in the accumulation area is too stiff, because flow parameters were derived from a calibration between glacier geometry and measured surface velocities in the ablation zone. Alternatively, the reconstructed ice thickness in the accumulation zone, inferred from the concept of plastic flow with a constant yield stress, may be an underestimate. Quite possibly the real ice thickness in the accumulation zone and the real ice volume in 2001 are much closer to the modelled estimate of 1.41 km³.

Fig. 6. Ice thickness distribution in 2001, derived from field measurements (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013) (a) and from the transient model run (b). The contours represent the difference in ice thickness (m) between the two methods.

A final comparison can be made between the surface mass balance predicted by the energy-balance model and the actual stake readings from our observations. RMSEs for individual balance years vary between 0.5 and 1.0 m w.e. a⁻¹ and are similar to those obtained during the model calibration by construction (Reference Nemec, Huybrechts, Rybak and OerlemansNemec and others, 2009). The largest discrepancies are found in the lower part of the glacier for the warmest years. However, the deviations have a stochastic component and did not prevent the model from correctly reproducing the recent glacier evolution as discussed above.

4.1. Climatic-forcing scenarios

With the past evolution of the glacier well constrained, the simulations can be extended into the future. To force the future mass balance we use climatic scenarios inspired by the CH2011 (2011) projections for the 21st century over Switzerland. These are based on output from global climate models (GCMs) that are coupled to European-scale regional climate models (RCMs). CH2011 considers three Inter-governmental Panel on Climate Change (IPCC) scenarios, of which two are non-intervention scenarios (A2 and A1B) and a third is a climate-stabilization scenario (RCP3D, emissions cut by 50% in 2050 compared with 1990). With respect to 1980–2009, the different GCM/RCM combinations and scenarios result in a mean temperature increase in 2070–99 in the range of +0.8 to +2°C (RCP3PD scenario) and 2.7°C to +5.0°C (A2 scenario). For the Alpine region, the warming in all RCMs is most pronounced in the summer months. Also the change in precipitation shows a seasonal character, with wetter winters (between +0% and +30% precipitation) and drier summers (–30% to 0%). This is in line with the results of the PRUDENCE project (Prediction of Regional scenarios and Uncertainties for Defining European Climate change risks and Effects; Reference Christensen and ChristensenChristensen and Christensen, 2007). Based on these climate projections, we construct nine schematic warming scenarios that broadly cover the range of CH2011 scenarios. We opted for schematic scenarios as the interest is in the main future climatic trends and not in a specific scenario-model set-up. As a control we also consider three scenarios with no temperature increase. The nine warming scenarios consist of a combination of three linear warming trends (+1.0°C, +2.5°C and +4.0°C of mean annual warming between 2010 and 2100) and three precipitation trends (neutral, wet and dry). The changes are applied with respect to the 2001–10 averages of the monthly mean temperature and precipitation data. For all scenarios, the summer (June–July–August) warming lies 25% of the annually averaged warming above the warming of the other seasons: for example, the +4.0°C scenario corresponds to a temperature increase of +4.75°C for the summer and +3.75°C for the other seasons. The wet scenario accounts for a linear 30% increase of winter precipitation by 2100, an increase of 15% for spring and autumn and no change in summer precipitation. In the dry scenario the summer precipitation linearly decreases by 30%, the spring and autumn by 15% and the winter precipitation remains unchanged. For the neutral precipitation scenario, future precipitation rates are kept constant at their 2001–10 level.

4.2. Projections of ice volume and ice extent

Figure 7 and Table 1 summarize the projected evolution of glacier volume. Under all scenarios a strong mass loss is expected. This includes the no-warming scenario, in which the glacier still loses between 24.4% and 33.4% of its volume by 2100 depending on whether precipitation rates increase, decrease or remain stable. This provides a clear manifestation of the current imbalance of Morteratsch glacier and of the time needed to equilibrate towards 2001– 10 mass-balance conditions. For the highest warming scenario considered here (+4.0°C), ice volume by 2100 is reduced by 80% relative to 2010. For all warming scenarios, ice loss continues beyond the 21st century. On the other hand, up to 2020 the glacier evolution does not depend on any particular scenario and is almost entirely conditioned by past climatic changes, in particular the post-1970 warming. Strikingly the effect of precipitation changes is small compared with the effect of temperature changes. Even the adopted increase in precipitation of the wet scenario is not sufficient to counterbalance further glacier retreat for the no-warming case. This relates to the fact that the annual mass balance of Alpine glaciers is mainly driven by their summer balance (Reference Vincent, Kappenberger, Valla, Bauder, Funk and MeurVincent and others, 2004; Reference Zekollari and HuybrechtsZekollari and Huybrechts, 2014).

Fig. 7. Future evolution of ice volume under different climate scenarios. The black line represents the reconstructed evolution over the past 50 years. The line colour denotes the temperature scenario, while the line style indicates the precipitation scenario. For a particular colour, the thick solid line represents the ‘neutral’ precipitation scenario, the thin dashed line represents the ‘wet’ scenario and the thin solid line the ‘dry’ scenario.

The corresponding glacier geometries are shown in Figure 8. By 2030 the glacier has retreated by ~1 km compared with 2010, and Vadret Pers is no longer a tributary of Vadret da Morteratsch. The separation of the two glaciers occurs in the model before 2020. The imminent disconnection of the two glaciers is clearly supported by field observations since 2000 and may perhaps take place as early as 2015 if a remaining frontal ice patch of unknown thickness connecting both glaciers melts completely. For both glacier tongues, the frontal thinning and retreat is almost entirely driven by the local surface mass balance as there is virtually no mass supply from upstream as flow velocities in the tongue area approach zero (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). Typical ablation rates at the tongue of Morteratsch glacier are 7–8 m w.e. a⁻¹ for the period 2001– 10, and so are observed ice-thinning rates. For the highest-warming scenario (+4.0°C by 2100), the mean ELA has risen by 40–80 m by 2030 to ~3100 m for Vadret da Morteratsch and 3120 m for Vadret Pers.

Fig. 8. Simulated glacier extent and ice thickness for different temperature scenarios (vs 2010) and time periods, assuming neutral precipitation. The ice-thickness scale is the same as that used in Figure 6. The thick black line delineates the DEM glacier outline in 2001, while the thin lines delineate 50 m ice-thickness intervals.

After 2050, the effect from considering different scenarios becomes evident. By 2065, the retreat is stronger for Vadret Pers than for Vadret da Morteratsch, mostly due to thinner ice, a higher ELA and a more negative mass balance in the frontal area. In the +4.0°C warming scenario, the ELA for Vadret Pers rises to 3330 m by 2065, while for Vadret da Morteratsch it rises to ~3315 m. Towards the end of the century the ELA rises to 3520 m for Vadret da Morteratsch and 3500 m for Vadret Pers. This reflects a similar exposure and orientation of the respective glacier sections in this elevation band. Glacier retreat continues up to the end of the 21st century. By 2100, there is a clear distinction between the different scenarios considered. For the two highest-warming scenarios, only ice at high elevations remains and the glacier complex starts to disintegrate into disconnected ice patches. Despite its slightly more negative mass balance, the tongue of Vadret da Morteratsch extends further north than that of Vadret Pers because of its higher ice thickness which takes longer to completely waste away.

Three-dimensional views of the various model simulations discussed above are displayed in Figure 9. In the scenarios with the strongest retreat, it is likely that proglacial lakes will form in bedrock depressions (Reference Zekollari, Huybrechts, Fürst, Rybak and EisenZekollari and others, 2013). Their effect is not taken into account here, but can potentially accelerate the glacier retreat by hydrostatic lifting of the glacier front. The formation of future proglacial lakes is also important from a natural-hazards perspective, as some of them form behind natural moraine dams that can easily break, potentially causing serious damage from flash flooding downstream (Reference Haeberli, Alean, Müller and FunkHaeberli and others, 1989; Reference Frey, Haeberli, Linsbauer, Huggel and PaulFrey and others, 2010). Judging from the glacier bed exposed so far, the lake that may form in the central Morteratsch area is, however, probably dammed by hard rock rather than by soft sediment, and so may not be subject to potential failure.

Fig. 9. 3-D view of the glacier from the north-northwest for different settings, scenarios and times. The domain ranges from 789 500 m to 796 000 m in the x –direction and from 137 500 m to 148 000 m in the y –direction following the CH1903 coordinate system. Blue areas denote local bedrock depressions that may hold proglacial lakes after the glacier has retreated.

4.3. Climate–geometry imbalance

Further insight into the future time-dependent response and glacier imbalance is provided by the time evolution of the mean specific mass balance. Figure 10 displays the spatially averaged annual surface mass balance for all future warming experiments with constant (2001–10 average) precipitation. The mean specific mass balance is a measure for the annual rate of mass loss and is a complex function of both the climate state and the actual glacier hypsometry (area– elevation distribution). For constant 2001–10 climate forcing (no-warming scenario), the expected quasi-exponential decay of the glacier imbalance is confirmed. It is caused by geometric adjustment that mainly entails loss of an excess glacier tongue that can no longer be supplied with sufficient ice mass from a reduced accumulation area.

Fig. 10. Future evolution of mean specific mass balance for the different climate scenarios, assuming average 2001–10 precipitation.

By the end of the 21st century, the glacier is still out of equilibrium with the 2001–10 climate, although an asymptotic evolution towards a zero mean specific mass balance is evident. The eventual steady state corresponds to a volume loss with respect to 2010 of between 25% and 35%, depending on the precipitation scenario, and a frontal retreat of almost 2 km for Vadret da Morteratsch and 1.5 km for Vadret Pers (see also Figs 7–9). In this experiment, the associated AAR steadily increases from 44% for the period 2001–10 to 57% in 2065 and 58% in 2100.

From the no-warming case, one can also derive a first estimate of the volume response time of Morteratsch glacier, as the time period during which a fraction of (1 – e⁻¹) or 63% of the imbalance in 2010 of –0.82 m w.e. a⁻¹, derived from the 2001–10 mean climatology, has been overcome. This timescale corresponds to 32 years but is not considered a clean measure of the response time of the current glacier as it also contains a contribution from the cumulated imbalance prior to the 2001–10 period.

For the other warming scenarios, the evolution of the mean specific mass balance results from an interplay between geometric adjustment and a continuously decreasing surface mass balance. In the +4°C scenario the effect of further warming initially dominates and the imbalance increases until 2052. After that period, the loss of surface area is more pronounced and this increases the mean specific mass balance faster than the warming decreases it. The complex evolution of the imbalance in a changing climate illustrates the difficulty of deducing a climatic signal from changes in glacier geometry. For this purpose the reference surface mass balance, the mass balance on a fixed surface at a certain reference date (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001; Reference Huss, Hock, Bauder and FunkHuss and others, 2012), should be used. Over the 21st century, the glacier’s reference surface mass in our experiments is found to increase continuously for all warming scenarios (i.e. it becomes less negative).

5.1. Comparison with other studies

Comparable studies simulating the future 3-D geometry of a mountain glacier have been performed by other researchers. Using a distributed accumulation and temperature-index melt model driven by RCMs to estimate ice-thickness changes, Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) attempted to project the future evolution of the Morteratsch glacier complex. In their modelling, similar climate scenarios were used but ice dynamics was not considered. Instead, they imposed a distributed ice-thickness change based on the observation that elevation changes of retreating glaciers are normally largest near the glacier terminus. Since Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) also carried out a model calibration over the 20th century, there is reasonable agreement between their projections and ours concerning the evolution of glacier length, area and volume until 2030. Beyond 2030, however, Reference Huss, Usselmann, Farinotti and BauderHuss and others (2010b) produce a more fragmented glacier retreat with more isolated patches of ice in the deeper bedrock depressions, which are largely absent from our work as analysed further below.

The relative volume loss of the Morteratsch glacier complex over the 21st century is generally smaller than that of other large Swiss glaciers (Reference Huss, Usselmann, Farinotti and BauderHuss and others, 2010b; Reference Farinotti, Usselmann, Huss, Bauder and FunkFarinotti and others, 2012; Reference Salzmann, Machguth and LinsbauerSalzmann and others, 2012). This difference depends on several factors, such as the sensitivity of the mass balance and differences in climatic scenarios, but is mainly dictated by the specific hypsometry of the glacier in question. Morteratsch glacier stands out in having a comparably large areal fraction in the elevation band between 3500 and 4000 m, which will remain in the accumulation zone even for a warming of +4°C. In another 3-D ice-dynamic simulation, Reference Jouvet, Huss, Blatter, Picasso and RappazJouvet and others (2009) find that Rhonegletscher almost entirely disappears by 2100 in their ‘median’ scenario (in between our +2.5°C and +4.0°C scenarios), but that is mainly because Rhonegletscher has its highest point at ~3600 m, eliminating an accumulation area of sufficient size. In another study with the same ice-dynamic model, Reference Jouvet, Huss, Funk and BlatterJouvet and others (2011) reported dramatic ice losses for Grosser Aletschgletscher by 2100. In their ‘ENSmed’ scenario, close to our +4.0°C scenario, <10% of the 1999 ice volume remains by 2100, less than the ~15% remaining for Morteratsch glacier. For Glacier de Saint-Sorlin (2700–3400 m) Reference Le Meur, Gerbaux, Schäfer and VincentLe Meur and others (2007) modelled a total disappearance of the glacier by the end of the 21st century under the IPCC’s SRESB1 scenario. Again, these differences reflect different area–elevation distributions. A common feature of all these studies, however, is the near disappearance of even the largest Alpine glaciers by the end of the 21st century, albeit for a higher-than-average warming scenario.

5.2. The role of ice dynamics

To further investigate the role of ice dynamics in our projections, we considered two additional experiments. These test the role of ice dynamics in different set-ups following the work of Reference Huybrechts and WoldeHuybrechts and De Wolde (1999). The dynamic response is the one discussed so far in which there is a full coupling between the mass balance and the ice flow. In the static response, the ice flow is kept constant at the year 2010 state (constant flux divergence in Eqn (6)) but changes in ice thickness can still feed back on the mass balance through the mass-balance elevation gradient. For the fixed-geometry response, the ice flow is also kept constant, and additionally changes in glacier geometry cannot feed back on the mass balance (through the mass-balance elevation gradient). Here only perturbations with respect to the 2001–10 mass balance, which are solely climate-driven, are integrated forward in time. This set-up enables us to distinguish between the role of ice dynamics on the one hand and the effect of the height–mass-balance feedback on the other hand.

Figure 11 shows the results. For both the no-warming and the +4°C scenarios the differences in total volume loss by 2100 are within 10% of each other for all three cases. This clearly shows that the dominant response is from changes in the surface mass balance, not from ice dynamics. The static and fixed-geometry cases are most closely aligned, indicating that due to the limited ice thicknesses the height–mass-balance feedback does not play an important role in the future response of an alpine mountain glacier. This is different from the response of ice sheets (Reference Huybrechts and WoldeHuybrechts and De Wolde, 1999), which have a much larger ice thickness.

Fig. 11. Future evolution of glacier volume under different model set-ups addressing the role of ice dynamics. The upper curves are for the no-warming scenario and the lower curves for the +4.0°C warming scenario.

As time evolves, the inclusion of ice dynamics is found to speed up the total volume loss, and this effect is much more pronounced in the +4°C scenario than in the no-warming scenario. Moreover, keeping the flux-divergence term constant in the static response has an opposite effect on the accumulation and ablation zones. This is demonstrated in Figure 12. Including ice dynamics increases the mass loss in the ablation zone and enhances the retreat of the glacier tongue. This is because the flux-divergence term can adjust to the thinning margin (it becomes less positive) so that less mass is transported to the same position in the tongue area. This is not the case for the static experiment in which the ice-transport term is fixed. In the accumulation zone, on the other hand, excluding ice-dynamic adjustment (static response) increases mass loss and enhances local thinning, so that the glacier disintegrates more easily into separate ice patches than in the dynamic response.

Fig. 12. Glacier geometry in 2100 in the +4.0°C scenario in an experiment including ice-dynamic adjustment (dynamic response, left panel) and excluding ice-dynamic adjustment (static response, right panel). The view is from the north-northwest and coordinates are given in the CH1903 system.