Dynamic modelling of future glacier changes: mass-balance/elevation feedback in projections for the Vestfonna ice cap, Nordaustlandet, Svalbard

Abstract

Future projections of the evolution of ice caps as well as ice sheets and consequent sea-level rise face several methodological challenges, one being the two-way coupling between ice flow and mass-balance models. Full two-way coupling between mass-balance models – or, in a wider scope, climate models – and ice flow models has rarely been implemented due to substantial technical challenges. Here we examine some coupling effects for the Vestfonna ice cap, Nordaustlandet, Svalbard, by analysing the impacts of different coupling intervals on mass-balance and sea-level rise projections. By comparing coupled to traditionally deployed uncoupled strategies, we prove that neglecting the topographic feedbacks in the coupling leads to underestimations of 10–20% in sea-level rise projections on century timescales in our model. As imposed climate scenarios increasingly change mass balance, uncertainties in the unknown evolution of the fast-flowing outlet glaciers decrease in importance due to their deceleration and reduced mass flux as they thin and retreat from the coast. Parameterizing mass-balance adjustment for changes in topography using lapse rates as a cost-effective alternative to full coupling produces satisfactory results for modest climate change scenarios. We introduce a method to estimate the error of the presented partially coupled model with respect to as yet unperformed two-way fully coupled results.

1. Introduction

The two main modes of response of land-based ice bodies to climate perturbations are changes in ice dynamics and climatic mass balance (e.g. Moore and others, 2013). Ice dynamics are characterized by the type of flow regime and the resulting discharge into the ocean, in case of tidewater glaciers. Climatic mass balance (CMB) is the sum of surface mass balance (primarily the difference between snow accumulation and melt and sublimation of snow and ice, i.e. ablation) and the internal mass balance, excluding mass changes through ice flow, but comprising melting and refreezing inside the ice body (Cogley and others, 2011). Both ice dynamics and CMB are coupled processes with interactions becoming more important on long timescales when, for instance, changes in ice-sheet topography may significantly affect both CMB and ice dynamics (Church and others, 2013). The most rapid changes so far observed in terms of impacts on sea level are due to rapid acceleration of fast-flowing outlets (e.g. Rignot and others, 2011; Shepherd and others, 2012; Moore and others, 2013). Understanding the CMB/ice-dynamics combined response to changing climate is crucial for predictions of future sea-level rise (Moore and others, 2011; Rignot and others, 2011; Dunse and others, 2012; Shepherd and others, 2012).

Various issues make reliable future predictions difficult, for example the understanding of the basal processes leading to observed accelerations or surges of fast-flowing outlets (Schäfer and others, 2014, and references therein), and the lack of a physics-based calving law properly formulated for three-dimensional (3-D) configurations (Benn and others, 2007; Åström and others, 2014; Cook and others, 2014). Another limiting factor is the correct implementation of the retreat of marine-terminating glaciers in some ice flow models (Gagliardini and others, 2013).

Here we focus on another sparsely addressed issue: coupling of realistic dynamic ice-sheet models to predictive climate models (Bindschadler and others, 2013; Moore and others, 2013). We focus on the feedback of the CMB to time-evolving ice-sheet topography. Until recently, climate models, which are expensive in CPU time, have usually been deployed independently of ice-sheet models. The CMB fields obtained are then used as climatic forcing for separate runs of ice-sheet models which simulate the dynamic response of the ice bodies. Such an uncoupled configuration is not able to account for feedbacks to the CMB expected as a response to the evolution of the surface topography (Church and others, 2013).

CMB may be calculated with physically based, sophisticated representations of the mass and energy budgets associated with snow and ice surfaces in regional or global climate models, such as RACMO (Van Pelt and others, 2012) or WRF (Skamarock and others, 2008); these are, however, uncoupled from ice flow models and thus based on fixed surface topographies. Alternatively, some very simplistic CMB ‘models’, based on an empirically derived positive degree-day (PDD) scheme, may be directly incorporated into the ice flow models, taking the time-evolving surface topography into account. In these cases the mass balance is parameterized as a simple function of precipitation and air temperature (Reeh, 1991; Edwards and others, 2014a).

Some coupled simulations have been conducted, mostly for the Greenland ice sheet, often using ice-sheet models with simplified physics and relatively coarse-resolution regional climate models. All share the conclusion that predictions using stand-alone ice-sheet models omitting climate/ice-sheet feedbacks, should be treated with caution, at least when large changes in the ice sheet occur. Ridley and others (2005) couple the Greenland ice sheet model (GISM) to a global circulation model and investigate a locally induced climate change. Driesschaert and others (2007) use a 3-D Earth System Model of intermediate complexity and show that climate feedbacks considerably enhance the greenhouse-gas-induced warming over Greenland. Mikolajewicz and others (2007) couple ice-sheet (SICOPOLIS; Greve, 1995), atmosphere–ocean general circulation and vegetation models and find pronounced effects on estimates of the future development of the Greenland ice sheet, but only small impacts on large-scale climate. Swingedouw and others (2008) found enhanced ice loss from Greenland when coupling models, with changes in ice topography and meltwater flux influencing ocean and atmospheric circulations as well as sea-ice distribution. Vizcaino and others (2010) found reduced shrinkage of the Greenland ice sheet with coupling, mainly caused by the effect of topographic changes on surface temperature. Goelzer and others (2013) use a higher-order, 3-D thermomechanical ice flow model and compare the use of a high-resolution regional climate model with a PDD formulation for the Greenland ice sheet. Interaction between surface mass balance and ice discharge limits the importance of outlet glacier dynamics with increasing atmospheric forcing, and the runs using the regional climate model produce a significantly higher sea-level contribution compared with using the PDD formulation. Following Helsen and others (2012), Edwards and others (2014b) use a simplified ‘coupling’ based on a set of gradients that relate surface mass-balance changes to height changes instead of full coupling and explore the modelled response of the Greenland ice sheet within various ice-sheet models observing an additional sea-level contribution due to the mass-balance/ elevation feedback. Helsen and others (2013) similarly use a two-way asynchronously coupled regional climate/ice-sheet model including physically realistic feedbacks between the changing ice-sheet topography and climate forcing to estimate the contribution from the Greenland ice sheet to the Eemian sea-level highstand – an example for an application when full coupling is especially challenging. In their ‘asynchronous two-way coupling’ they update the surface topography in the climate model every 1.5 ka, and in between apply a correction using the gradient relating surface mass-balance changes to height changes.

This study aims to provide insights into the importance of full two-way coupling between a CMB model and a full-Stokes ice-sheet model run on an ice cap. Since full coupling is technically complex and computationally expensive, we focus on Vestfonna, a medium-size ice cap covering ∼2340 km² (Braun and others, 2011) on Nordaustlandet, the second largest island of the Svalbard archipelago. To assess the effect of two-way coupling compared to uncoupled simulations, we compare the surface elevation and ice volume change of the ice cap as well as integrated CMB changes in response to various century-long scenarios of climate forcing utilizing different coupling intervals (including no coupling) and a lapse rate approach.

The topographic feedback affects both ablation and accumulation via air temperature, precipitation and incoming radiation. Melting of the ice lowers the surface elevation and thus exposes the ice cap to even higher air temperatures. This, in turn, leads to changes in the distribution of glacier facies and a general lowering of the surface albedo in the affected areas and thus higher amounts of absorbed global radiation. Hence, under the influence of a warming climate, both processes create a self-sustaining, positive feedback loop regarding ablation. Precipitation tends to increase with surface elevation, and thus a negative feedback on accumulation occurs for lower ice-cap topographies. As refreezing is affected by both ablation and accumulation, substantial feedbacks can also be expected for this mass- and energy-balance-relevant process but quantitative predictions are not feasible without considering dedicated model calculations.

By lapse rate approach we mean forcing of the ice flow model by uncoupled CMB input adjusted approximately for elevation changes using an elevation lapse rate similar to the gradient approach of Edwards and others (2014b). Such a lapse rate approach could be seen as a simple way of taking into account topographic feedbacks between the ice flow model and the CMB model which allows the climate or even the CMB models to be run on larger grids with coarser resolution uncoupled from sophisticated ice flow models that run at finer resolution.

We use output from the MIROC-ESM Earth System Model representing the four representative concentration pathways (RCP) 2.6, 4.5, 6.0 and 8.5 (Moss and others, 2010) created as part of the Coupled Model Intercomparison Project 5 (CMIP5) as climate forcing for a CMB model that has been designed to cope with local characteristics at Vestfonna (Möller and others, 2013; Möller and Schneider, 2015). The ice flow model used in our study is Elmer/Ice (Gagliardini and others, 2013). The dynamics of the fast-flow glaciers are parameterized by results from inverse modelling (Schäfer and others, 2014) using satellite-derived velocity fields.

We first outline the data and observations for Vestfonna. The flow model as well as the climate scenarios and the CMB model together with the applied coupling schemes are briefly summarized at the beginning of Section 3, followed by the different simulations we performed. We then suggest some alternatives when full coupling is not possible.

2. Location and Data

2.1. Research location

The European Arctic archipelago of Svalbard has ∼34 560 km² of glaciated area (Moholdt and others, 2010a). Nordaustlandet, the second largest island of the archipelago, is, to a large extent, covered by the two major ice caps Austfonna (8000 km²; Moholdt and others, 2010b) and Vestfonna (2340 km²; Braun and others, 2011; Fig. 1a). In this study we focus on the ice cap Vestfonna (VSF), which covers the western part of Nordaustlandet. Its dome-shaped ice body reaches 630 m at its highest point (Braun and others, 2011), Ahlmann summit, which is located in the centre of an extensive east–west-stretching ridge. From a minor summit in the eastern part of the ice cap a second ridge extends northwards. Apart from these main terrain features, the ice cap is dominated by large outlet glacier basins that drain into the surrounding fjords. Subglacial bedrock lies below sea level around some of the periphery (elevations down to −160 m a.s.l.), while in the central parts bedrock elevation increases to ∼410 m a.s.l. With a mean ice thickness of 186 m the ice cap holds a volume of 442 km³, which is relatively small for its surface area (Petterson and others, 2011).

Fig. 1. (a) Geographical location of Vestfonna. The red dots mark the location of MIROC-ESM gridpoints; the bigger dot is the location used for the cloud data. Image courtesy of NASA. (b) Initial surface elevation and (c) climatic mass balance in 2006/07 at the start of the simulations. The contours in (b) show the ice thickness. The grey cross in (b, c) marks the automatic weather station (AWS) in the west of Vestfonna, used for downscaling air-temperature and precipitation data. The outlet glacier Franklinbreen and the cross section used in Figure 11 are shown in (c).

Most calving fronts have experienced continuous retreat over the past three decades without significant change in flow velocities or any clear periodic, seasonal accelerations (Braun and others, 2011; Pohjola and others, 2011). An exception to this general behaviour is the largest outlet glacier, Franklinbreen, which drains northwestwards into Lady Franklinfjorden. This glacier has undergone a small but continuous advance combined with substantial speed-up (Braun and others, 2011; Pohjola and others, 2011), suggesting a radical change in flow-velocity regime but with no clear indication of a typical surge (Schäfer and others, 2014). However, past surges have been reported for Franklinbreen and other outlets of Vestfonna (Błaszczyk and others, 2009).

A mass-balance year at VSF typically lasts from the beginning of September until the end of August, with the accumulation period covering the first 9 months of this period (Möller and others, 2011a). The ice cap features a large accumulation area in its central part and a surrounding ablation area below the equilibrium line between 300 and 450 m a.s.l. depending on annual mass balance (Möller and others, 2013). Accumulation shows high temporal variability, with changes of up to 100% between individual years (Beaudon and others, 2011). Spatial variability across the ice cap is governed by terrain elevation and snowdrift (Möller and others, 2011b; Sauter and others, 2013).

Numerical modelling studies by Möller and others (2011a, 2013) suggest that the CMB of Vestfonna has been in a slightly positive state over recent decades. For mass-balance years 1979/80 to 2010/11 a mean CMB rate of +0.09 ± 0.15 m w.e. a⁻¹ was obtained, with annual balances showing a slightly negative but insignificant trend that led to an increase in equilibrium-line altitude over time (Möller and others, 2013). This increase is likewise not statistically significant. Following this development, the CMB dropped to a mean rate of −0.02 ± 0.20 m w.e. a⁻¹ over the period 2000/01 to 2008/09 (Möller and others, 2013). This is consistent with remote-sensing-derived geodetic mass balance of the ice cap (Moholdt and others, 2010b; Nuth and others, 2010).

2.2. Data basis and preparation

The study requires topographic data of VSF, i.e. bedrock and surface elevations. Bedrock elevation (Petterson and others, 2011) is kept constant while surface topography is updated according to the modelling experiments. As initial surface topography we use the digital elevation model (DEM) from the Norwegian Polar Institute (NPI) (1 : 100 000, 1990, UTM zone 33N, WGS 1984), which is based on topographic maps derived from aerial photography, completed with the International Bathymetric Chart of the Arctic Ocean (Jakobsson and others, 2008) for surrounding sea-floor. Details of the surface velocity data used for inferring the basal friction parameters are provided by Schäfer and others (2014).

To run the CMB model, daily air temperature, precipitation and cloud-cover data from the period September 2006 to August 2100 representing RCPs 2.6, 4.5, 6.0 and 8.5 (Moss and others, 2010) are used (Fig. 2). These data are taken from the Earth System Model MIROC-ESM (Watanabe and others, 2011). MIROC-ESM provides data at a resolution of 2.8125° × 2.8125°. According to a comparison of ten different ESM and global climate models carried out by Möller and Schneider (2015), MIROC-ESM predicts the highest air temperature increase in the study area over the 21st century. It was selected as forcing in this study as it leads to the strongest changes of the ice masses of Vestfonna by 2100 and thus potentially also produces the largest differences among our coupling experiments.

Fig. 2. Downscaled annual values of (a) mean annual air temperatures, (b) annual precipitation sums and (c) annual mean cloud cover on data from four gridpoints (Fig. 1a) from the MIROC-ESM Earth System Model under RCP 2.6, 4.5, 6.0 and 8.5. These are used as input to the climatic mass-balance model.

For the simulations in this study, we use daily data of air temperature and precipitation from the four ESM gridpoints closest to VSF (Fig. 1a) that were downscaled by Möller and Schneider (2015) to fit local conditions on the northwestern slope of VSF as represented by adjusted European Centre for Medium-Range Weather Forecasts ERA-Interim reanalysis data described in Möller and others (2013). Air temperature downscaling was based on a combination of multiple linear regression and variance-inflation techniques (Karl and others, 1990; Huth, 1999; von Storch, 1999; Möller and others, 2013). Precipitation downscaling combined multiple linear regression and local scaling techniques (Widmann and others, 2003; Möller and Schneider, 2008).

The air temperature downscaling for data of RCP 2.6 (4.5, 6.0, 8.5) resulted in a monthly root-mean-square error (RMSE) of 2.42 K (2.36, 2.23, 2.29 K) for months of the accumulation season, i.e. September to May, and an RMSE of 0.54 K (0.62, 0.50, 0.64 K) for months of the ablation season, i.e. June to August. The downscaled precipitation data of RCP 2.6 (4.5, 6.0, 8.5) show a monthly RMSE of 2.3 mm (2.4, 2.4, 2.3 mm) independent of mass-balance season.

The derived quantities, cumulative positive degree-days, cumulative snowfall and rain–snow ratio, are needed as additional input for the surface-albedo module. They are calculated from the downscaled air temperature and precipitation data according to Möller (2012). Daily cloud-cover data are directly taken from the gridpoint closest to VSF (79.53° N, 19.69° E); the other three gridpoints enclosing VSF are too far away to be relevant.

We use ERA-Interim based data (Möller and others, 2013) for mass-balance years 1996/97–2005/06 as the reference control for this study. These data also serve as the reference during GCM downscaling.

3. Methods

In our application, a CMB model links the climate model to an ice flow model. The models are run coupled, accounting for the topographic feedback. Note that, for larger ice bodies, the effect of changing ice-sheet topography should be fed back into the climate model and not only into the CMB model. Reasons for this include changes in the surface topography and surface albedo that could alter regional-scale air temperature, precipitation and atmospheric circulation; runoff and thus the freshwater input into the ocean that could affect ocean circulation (Ridley and others, 2005); and sea-surface temperature/sea-ice distributions (Swingedouw and others, 2008; Vizcaino and others, 2010). However, due to the relatively small size of VSF, we neglect this topographic feedback.

We first describe the ice flow model used to calculate the evolution of the VSF ice cap during the next century under different climatic scenarios. Then these scenarios are presented together with the applied climatic mass-balance model. Thereafter the implemented coupling of these two models is laid out. Finally, we describe the set-up of all our simulations before outlining some alternatives to full model coupling.

3.1. Model description

3.1.1. Ice flow model

The model equations are solved numerically with the Elmer/ Ice ice flow model, which is based on the open-source multi-physics package Elmer developed at the CSC – IT Center for Science, Espoo, Finland, and uses the finite-element method (Gagliardini and others, 2013). Here we use the methods described by Schäfer and others (2012, 2014).

The ice is modelled as a nonlinear viscous incompressible fluid flowing under gravity over a rigid bedrock. At the quasi-static equilibrium the force balance is described by the Stokes equations expressing the conservation of linear momentum and mass. Here we assume that ice is an isotropic material and hence its rheology is given by Glen’s law (Glen, 1952; Nye, 1952), linking deviatoric stresses and strain rates. The temperature dependency of the viscosity is described by an Arrhenius law (Schäfer and others, 2014).

The evolution of the free surface, S, is governed by the kinematic boundary condition

(1)

where b is the climatic mass balance (expressed in ice equivalent) and v_{x ,y,z} are the components of the velocity vector, .

S is approximated as a stress-free surface. On the lateral boundaries, the normal stress is set to the water pressure exerted by the ocean for elevations below sea level, otherwise we prescribe a stress-free condition.

In the absence of a calving law, we chose to prescribe fixed margins at marine outlets. Hence the lateral mass loss is a consequence of the flow solution rather than a property prescribed by a physical calving law. On VSF, lateral mass loss accounts for roughly one-fifth (or less) of the total volume loss. Figure 3 shows the time evolution of the lateral mass loss and the contribution of CMB for a fully uncoupled simulation under RCP 2.6. The lateral mass loss remains fairly constant over time and is dominated by increasing loss through CMB. Retreat of land-terminating margins is modelled by thinning of the ice in combination with a minimum flow depth of 10 m (smaller values are considered to be deglaciated).

Fig. 3. Time evolution of the lateral mass loss and loss introduced by CMB during a fully uncoupled simulation under RCP 2.6.

At the lower boundary, B, a linear friction law (Weertman law) in the form of a Robin boundary condition (Greve and Blatter, 2009) is imposed:

(2)

where and are normal and tangential unit vectors ( being in the tangential plane aligned with the basal shear stress), β is the basal friction parameter and and are the basal velocity and stress components parallel to the bed at the base. We assume zero basal melting .

The basal friction parameter field, β(x, y) is inferred in this study from surface velocities using an inverse method following the approach of Arthern and Gudmundsson (2010) as detailed in Schäfer and others (2014). At present it is impossible to define a time-evolving basal friction parameter field because of poorly understood physical processes. Hence we assume a temporally fixed basal friction parameter distribution, thereby reducing our investigation to a sensitivity study rather than a future projection. However, basal friction parameters based on data assimilation for two distinct time periods (1995 and 2008) are available, which introduces some temporal variability of this important parameter and allows for changes in one of the outlet systems (Franklinbreen). We mostly use the β distribution obtained for 1995, which in the area of Franklinbreen leads to much lower velocities than for 2008.

Our simulations are not thermomechanically coupled, because a correct, coupled spin-up is impossible to perform (Schäfer and others, 2014). This is a reasonable simplification to make for VSF given our aim of examining century-scale CMB feedback effects, and additionally saves computation time. Inside the ice body, the ice temperature is assumed to follow a linear depth-dependent temperature profile

(3)

where q _geo is the geothermal heat flux, assumed to be 40 mW m⁻², κ = 2.072 W K⁻¹ m⁻¹, a representative heat conductivity of ice for the ice temperature range, and is the ice depth (vertical distance to the surface at a given location in the ice body). We use the idealized case of surface temperatures adjusting instantly to the time-evolving air temperatures as prescribed in the climatic mass-balance model (detailed below). Additionally, we run some of our simulations keeping the ice temperature distribution fixed to its initial distribution at the beginning of the century. A short relaxation of the surface is made before proceeding with future simulations, allowing the mesh to evolve and to be internally consistent with respect to flow and ice temperature fields.

Mesh resolution is a critical factor in the simulation. Following Schäfer and others (2014) we use an anisotropic mesh utilizing the fully automatic, adaptive, isotropic surface re-meshing procedure YAMS (Frey, 2001) for the ice flow model. An initially homogeneously spaced two-dimensional footprint mesh was established using Gmsh (Geuzaine and Remacle, 2009) following the glacier outline of the NPI DEM, and refined with a metric based on the Hessian matrix of the observed 1995 surface velocities (Pohjola and others, 2011). Horizontal resolution varies between 250 and 2500 m. The resulting footprint mesh is vertically extruded to represent the 3-D shape of the ice cap using ten equidistant terrain-following layers.

3.1.2. Climatic mass-balance model

The CMB, b, which is surface mass balance plus refreezing (Cogley and others, 2011), is modelled using the temperature–net-radiation index model of Möller and others (2013). This model builds on that developed by Möller and others (2011a) but incorporates monthly updated surface albedo-field calculations according to Möller (2012). Calculations of accumulation, ablation and refreezing are done daily while the surface-albedo fields are updated on a monthly basis.

Surface accumulation, c(z), is calculated as an altitude profile according to

(4)

with P _ESM being the downscaled, daily MIROC-ESM precipitation at sea level and T(z) the downscaled, daily MIROC-ESM 2 m air temperature (Fig. 2) that is distributed over terrain elevation using a constant linear lapse rate of −7.0 K km⁻¹ (Möller and others, 2011a). The empirical parameters f _P1 (0.99 ×10⁻⁵ m⁻²) and f _P2 (0.79 × 10⁻² m⁻¹) were calibrated by Möller and others (2011a) using in situ snow-pit data from Möller and others (2011b). The transition from liquid to solid precipitation between +2 and 0°C is realized using a hyperbolic function (Möller and others, 2007).

Surface ablation, a(x, y, z), is calculated for days with positive mean air temperatures as a spatially distributed quantity according to

(5)

with α(z) being the modelled surface-albedo profile and R(x, y, z) the modelled global radiation. The empirical parameters, f _T (1.736 mm w.e. K⁻¹ d⁻¹) and f _R (0.14 mm w.e. W⁻¹ m² d⁻¹), were calibrated by Möller and others (2013) using in situ point mass-balance data from repeated stake measurements.

Surface albedo, α(z), is calculated as an altitudinal profile according to

(6)

with Θ(z) being a sigmoid function of rain–snow ratio and Ψ(z) a linear function of cumulative snowfall and cumulative PDDs. Further details of the albedo calculation are provided by Möller (2012).

Global radiation R(x, y, z) is calculated as a spatially distributed quantity according to

(7)

with R _S(x, y) being cloud-cover-reduced direct solar radiation calculated following standard solar geometry rules (Möller and others, 2011a, and references therein). Multiple scattering and reflection, R _M(x, y, z), is calculated using an empirical approach that was developed and calibrated for application at VSF (Möller and others, 2011a). Refreezing, r(x, y, z), is realized by applying the P _max approach of Reeh (1991). Here we use P _max = 0.9 (Möller and others, 2013), which means that ablation is considered to refreeze within the glacier until 90% of the previous winter accumulation is reached. Ablation that occurs after this proportion is reached leaves the glacier as runoff. In contrast to the unstructured mesh for the ice flow simulations, the CMB model variables are represented on a regular grid with 250 m grid spacing.

3.1.3. Model coupling

The ice flow model is run with a time step of 1 year. The CMB model produces daily CMB fields from which annual sums for each mass-balance year (i.e. September to August) are derived. Both models are run uncoupled as well as with different coupling intervals: 50, 25 and 10 years. Full coupling (i.e. at each time step) could be done in principle, but was not implemented because of difficulties in porting the codes to a common platform. However, we present alternative methods to estimate possible fully coupled results from our coupling experiments.

The ice flow model is run over the entire model domain, which remains unchanged during the simulations. An ice mask delineates the glaciated area of the model domain (minimum flow depth 10 m), for which the CMB fields are calculated. Positive glacier-wide CMB that may lead to glacier advances on longer timescales are limited to a couple of mass-balance years at the beginning of the simulation period. Hence, no advance outside the area given by the ice mask needs to be modelled.

In the uncoupled run, the CMB model uses the initial surface topography and ice mask throughout the full 94 years (2006–2100) of simulation. After the CMB model has been run over the simulation period, all calculated annual CMB fields are passed to the ice flow model. The ice flow model then calculates evolving surface topographies and ice masks over the simulation period. In the coupled runs, the surface topography and the ice mask used by the CMB model are updated after each coupling interval with the values calculated by the ice flow model. Thereafter, the CMB fields for the following years are calculated and provided to the ice flow model.

For this exchange of variables between the two models an interpolation between the individual model grids has to be established. To that end we utilize an inverse distance weighting method (Shepard, 1968) using the closest gridpoints for distributing the CMB variable values required for the ice flow model over the unstructured mesh. Interpolation of values from the unstructured mesh back to the regular grid utilizes a method that first finds the corresponding element containing the coordinates and then uses the finite-element basis functions in order to interpolate the required value.

By feeding back the evolved surface topographies to the CMB model, the coupling process leads to considerable changes in the calculated CMB components: precipitation and thus accumulation; air temperature and surface albedo and hence ablation. Precipitation increases with terrain elevation according to a second-order polynomial function. This means that a lowering of the glacier surface results in decreased precipitation and thus accumulation. Because of the polynomial shape of the precipitation gradient these effects are larger in the upper parts of the ice cap than at its margins. Air temperature decreases linearly with terrain elevation according to a constant lapse rate. Glacier-surface lowering would thus result in higher air temperatures and thus increased ablation rates, especially if the lowering leads to a more frequent occurrence of melt days, i.e. days with above-zero mean air temperatures. Surface albedo, which depends on a complex function of precipitation- and air-temperature-related quantities, would also decrease with a lowering of the glacier surface. Global radiation is the only meteorological input variable that does not have a simple elevation dependency, and direct solar radiation is independent of elevation. The spatial variability of global radiation on a regional-scale surface depends mostly on local aspects and slopes. Since the overall dome-like shape of the ice cap is largely maintained over the simulation period, changes of aspect and slope are relatively small. Topographically induced changes in global radiation and the resulting ablation changes are thus of minor importance compared with other CMB-relevant input variables. Refreezing, which is determined by all meteorological input variables, shows considerable variability due to changing topography because of its complex relation to snow distribution and ablation. However, any real quantification is difficult because of the complex, nonlinear and temporally varying dependency of refreezing on climate variables.

3.2. Set-up of simulations

The simulations we conducted are summarized in Table 1. We compare couplings at 50, 25 and 10 year intervals with uncoupled simulations. In addition, some of the runs were carried out with different β fields. All simulations span the period 2006/07–2099/2100, start from the same initial surface elevation and CMB as shown in Figure 1b and c and were conducted for each of the RCPs 2.6, 4.5, 6.0 and 8.5. Unless indicated otherwise, the β field from 1995 and time-evolving surface temperatures are used.

Table 1. Overview of conducted simulations. The ‘Coupling_ID’ column indicates the coupling interval (c0 no coupling, c1 coupling every 50 years, c2 coupling every 25 years and c3 coupling every 10 years). The estimated full two-way coupling results (C4) are included for completeness. lr designates the simulations with a lapse rate. The ‘β’ column indicates which basal friction parameter field, and the last column which temperature approximation, is used

3.2.1. Control runs

Control runs (simulations labelled control) have been conducted to check for model drift arising from model imbalances. They are driven by pre-simulation conditions represented by the mean 1996/97–2005/06 CMB and surface temperature fields. The underlying pre-simulation climate is based on adjusted ERA-Interim data (see Section 2.2). The drift is very small, so no corrections were needed (see Section 4). The control runs were conducted with both β fields to gain a first impression of the importance of variations in this parameter (Fig. 4).

Fig. 4. Surface elevation at the end of the control runs (constant climate) with two different β fields (a, b), as well as the surface elevation difference between these two runs (c). As in Figure 1 the ice thickness is given by coloured contours.

3.2.2. Uncoupled simulations over 94 years

Completely uncoupled simulations (simulations c0) were conducted for all RCPs as reference runs. In all cases VSF retreats and thins: the stronger the changes in the applied climate forcing, the more pronounced the ice volume changes become (Fig. 5, upper row). This is also visible in the surface velocity distributions (Fig. 6). Velocities decrease drastically as the ice cap thins more under the stronger applied climate forcings. Since the lack of a time-evolving basal friction parameter law is an important factor for future predictions, we assess the impact by conducting the same runs with the 2008 β field. For RCP 8.5, where we expect the strongest changes, a different approximation for the ice temperature was used to assess the impact of our approximation for evolution of the ice temperature: ice temperatures were kept fixed at the initial distribution compared with ice temperatures that are instantly adjusted to evolving air temperatures. Results do not indicate large sensitivity and are not shown here.

Fig. 5. Final surface elevations at the end of the century for all RCPs and different coupling intervals. Each row represents one coupling interval (c0, c1, c2, c3), each column one RCP (2.6, 4.5, 6.0, 8.5).

Fig. 6. Surface velocities at the beginning of the simulation and at the end of the uncoupled simulations c0 for all RCPs. As in Figure 5, grey represents the deglaciated area and light blue the 50 m ice thickness contour. The scale is cut at 200 m a⁻¹. Velocities are decreasing as the ice cap thins under more negative CMB.

3.2.3. Simulations with intermediary couplings at 50, 25 and 10 year intervals

All RCP simulations were performed with coupling at intervals of 50 years (c1), 25 years (c2) and 10 years (c3). The simulations lead to substantially different CMB evolutions over the 21st century, with annual CMB being systematically more negative than in c0 (Fig. 7). Similar observations can be made for the glacier geometry (Fig. 5).

Fig. 7. Annual glacier-wide climatic mass balances for the simulations c0 (yellow), c1 (red), c2 (green) and c3 (blue) for each of the four RCPs. Coupling dates are indicated by colour-coded vertical dashed lines and labelled at the top.

3.3. Alternatives to full two-way coupling

The aim of fully, i.e. annually, two-way coupled simulations incorporating mass-balance models on the one side and ice dynamics models on the other is to present reliable estimates of glacier-wide CMB and, most important, to produce an optimal ice volume change estimate for the end of the modelling period and hence reliable sea-level rise (SLR) equivalents. However, a full two-way coupling is usually not feasible. Here we outline two alternatives when either no coupling or only long coupling intervals are used.

3.3.1. Using a lapse rate as an alternative to real coupling

A simple way of implementing an elevation feedback when using output from large climate models without the possibility of implementing coupling at all is a lapse rate approach (Edwards and others, 2014b). In this case uncoupled CMB fields are provided by the climate models (or here a separate CMB model), but are corrected with a lapse rate for the changes in topography when being fed into the ice-sheet model as climatic forcing.

The simplest possible approach is to calculate a mass-balance lapse rate, λ, that is constant over time, from the fully uncoupled CMB fields, , i.e.

(8)

where is the surface elevation change from the initial surface elevation at the location at time t. Here we use a value of λ(= 0.0047 m ice eq. a⁻¹ m¹) simply determined from the 2006/07 CMB, in order to provide an approximation to the effect of this approach compared with coupling. This is only a crude approximation. Plotting b(z) _{c 0} for different times and RCPs shows that a time-independent λ, or even a linear relationship b ∼ z, does not perfectly reflect the elevation dependence of the modelled CMB distributions (Fig. 8). Despite this limitation, simulations with all four RCPs (simulations lr) were conducted and produce, at least for the first half of the simulation period, results fairly close to those of the c3 simulations that represent our shortest coupling interval (Fig. 9). For more advanced simulation times the deviations from the shortest coupling interval runs and the lapse rate approach become larger with increasing climate scenarios, suggesting that improvements need to be considered.

Fig. 8. Elevation dependency for different CMB fields: CMB in 2006/07 at the beginning of the simulations, and the uncoupled CMB for 2099/2100 at the end of the simulations for all four RCPs. The lines represent the linear fits used to determine lapse rates.

Fig. 9. Time series of the volume loss for the different simulations in mm global sea-level rise equivalent. A control run (grey line running close to bottom of plot) as well as simulations with a different basal friction field (crosses) are also shown.

3.3.2. Estimation of full two-way coupling results from simulations with longer coupling intervals

From the calculated glacier-wide CMB values at the end of the modelling period it is evident that a strong relationship exists between the length of the coupling interval and the associated deviations in glacier-wide CMB (Fig. 10). Starting from this, we present a simple but effective method to predict the glacier-wide CMB at the end of the modelling period as it would be obtained by a fully two-way coupled modelling architecture. When considering this method it should be kept in mind that no fully two-way coupled experiment was carried out in this study, so the reliability of the proposed method cannot be evaluated in any quantitative way. However, a qualitative evaluation, i.e. visual inspection, of the results (Fig. 10) suggests that the proposed method leads to reasonable estimations of fully coupled results. Extending this approach it is further possible to derive ice-volume changes and thus SLR equivalents as they would result from a full two-way coupled modelling.

Fig. 10. Method of estimating fully two-way coupled glacier-wide CMB (top row) and sea-level rise SLR (bottom row) in 2099/2100 from the simulations with four coupling intervals (c0, c1, c2, c3) performed. Differences relative to uncoupled CMB or SLR simulations (%) are indicated by the black open squares. A second-order polynomial (black equation and line) is fitted to the c0, c1, c2 data points, and extrapolated through simulation c3 to the hypothetical C4 simulation (dashed line). The error range (grey wedge) is determined by the difference in the CMB or SLR between actual simulation c3, and its polynomial extrapolation (∊_{c 3}, grey equation). The final CMB or SLR deviations, including error range for the hypothetical C4 simulation, are shown in red.

As a prerequisite, we calculate the final, i.e. 2099/2100, deviations in glacier-wide CMB of the c1, c2 and c3 simulations from the uncoupled glacier-wide CMB of the c0 simulation. Afterwards, a five-step procedure is used to estimate the glacier-wide CMB in 2099/2100 as it would have been obtained by an annually coupled (C4) simulation (Fig. 10a).

• 1. The coupling intervals of the c0, c1 and c2 simulations (94, 50 and 25 years) are related to the respective 2099/2100 CMB deviations by fitting a second-order polynomial to the three data points (black formulas in Fig. 10).

• 2. The residual between the polynomial fit at a 10 year coupling length and the 2099/2100 CMB deviation of the c3 simulation (∊_{c 3}) is calculated.

• 3. ∊_{c 3} is assumed to determine an error range around the polynomial fit. The error range linearly increases from zero at c2 with decreasing length of the coupling interval through the residual at c3 to the 1 year coupling interval of the hypothetical C4 simulation.

• 4. The C4 2099/2100 CMB deviation is given by the value of the polynomial fit at a 1 year coupling interval and the associated uncertainty from step 3 above.

• 5. The C4 2099/2100 glacier-wide CMB value is finally calculated from the C4 CMB deviation and the un-coupled c0 glacier-wide CMB.

This procedure is carried out separately for all four RCPs and can likewise be transferred to ice-volume changes in terms of SLR equivalents (Fig. 10b). Figure 10 shows the estimation of fully two-way coupled results for both glacier-wide CMB and SLR in 2099/2100 for all four RCPs.

4. Results and Discussion

4.1. Control runs and configuration tests

A series of control runs (simulations control; constant climate) and uncoupled runs (simulations c0; with changing climate) were carried out to assess model sensitivities of the ice flow computation. Control runs with both β fields showed that model drift was very small (volume changes of ∼0.01 mm SLR equivalent compared with climate-change induced values roughly between 0.3 and 1 mm; Table 2; Fig. 9), and no corrections were applied. Surface topographies resulting from runs using the different β fields showed that local differences in surface elevation were similar to or even higher than the total surface elevation change observed over the whole simulation period (Fig. 4).

Table 2. Volume change, ΔV _s, during the century simulations in mm global sea-level rise. ‘Abs.’ and ‘Rel.’ are the ΔV _s differences with respect to c0 for each RCP. ΔCMB is glacier-wide, negative deviation relative to c0. Simulations with different basal friction parameter fields did not differ significantly from the corresponding c0. The C4 rows are the estimates for full two-way coupling, and the lr rows are estimates based on lapse rate. The last column indicates the ice-cap volume, ΔV _i, at the end of the simulation relative to the initial volume of 4.75 × 10¹¹ m³ ice

This picture is strengthened by modified c0 simulations using the two different β fields. Switching from one β field to the other leads to higher or lower surface elevations over most of the ice cap; the variations are most pronounced in the Franklinbreen area (Fig. 11), independent of the applied RCP. These surface elevation differences are similar in magnitude to those when feedback is neglected. However, since they are not homogeneous over the ice cap, their variability is at least partly compensated when the total volume loss of the ice cap is calculated, and lies between 0.02% and 0.1% of the uncoupled run with the usual β field (difference between crosses and yellow lines in Fig. 9). The stronger the climate change, the less relevant the absolute surface elevation changes caused by changing β fields become for projections into the future, which is similar to simulations for Greenland’s outlet glaciers (Goelzer and others, 2013).

Fig. 11. Cross sections through the ice cap (shown in Fig. 1c), with bedrock shaded brown showing the initial topography, and that at the end of simulations under different RCPs using different coupling intervals, lapse rate parameterization, and the alternative basal friction parameter β field from 2008. The two control run results are shown in the RCP 2.6 panel.

Simulations c0 and lr were run for RCP 8.5 using an intime constant surface-temperature distribution in addition to the usual surface-temperature distribution. Compared with continuously adjusting to air temperatures prescribed by the CMB model, a small and fairly constant shift of surface elevations occurs over the whole ice cap (not shown). This observation is independent of whether the runs are conducted uncoupled or using the lapse rate approach. Hence, our results with respect to the spatial variations of surface-elevation changes and the absolute changes in volume loss are insensitive to the surface temperature approximation. For RCP 8.5 we observe relative volume changes for c0 of 8.1% and for lr of 7.4% representing the adjusting and the constant surface-temperature approximations respectively. Simulations with adjusting surface temperatures slightly overestimate the volume deviations when neglecting the coupling feedback compared with those where surface temperature is kept fixed.

The real initial surface-temperature field differs anyway from our idealized surface-temperature profile. From Schäfer and others (2014) we conclude that over most of the ice cap the real englacial temperature is higher than our initial temperature field. Both friction heat and firn heating introduce more heat into the ice cap, which is then redistributed by advection. Hence, the simulations with constant ice temperatures certainly represent a simulation with a too cold ice body, while those with adjusting ice temperatures may represent a too warm ice body for at least the later portion of the simulation time. We expect the real ice temperatures and related volume changes to lie somewhere between our two temperature approximations.

4.2. Time series of glacier-wide CMB, ice thickness and SLR

The evolution of annual, glacier-wide CMB shows substantial differences between the various simulations. In simulation c0, the glacier-wide CMB drops to values between −1.5 m w.e. a⁻¹ for RCP 2.6 and −4.9 m w.e. a⁻¹ for RCP 8.5 (Fig. 7). At the end of the century, the lower parts of the ice cap thus experience CMB values of −7 m w.e. a⁻¹ for RCP 8.5, and only for RCP 2.6 does a small area remain in the centre of the ice cap with CMB values close to 0 (Fig. 12; decadal average at the end of the simulation period).

Fig. 12. Climatic mass balance at the end of the simulation (decadal average) for the different RCPs without model coupling, i.e. when neglecting topographic feedback. Note the different colour scale from Figure 1 because of the absence of an accumulation area at the end of the century.

The coupled simulations (c1 to c3), in contrast, lead to considerably different CMB evolutions over the 21st century, with annual CMB being systematically more negative than in simulation c0 (Fig. 7). The shorter the coupling intervals or the further into the simulation, the higher the CMB deviations become (Fig. 13). In simulation c3, the annual glacier-wide CMB drops to values between −1.7 m w.e. a⁻¹ for RCP 2.6 and −6.1 m w.e. a⁻¹ for RCP 8.5 (Fig. 7). At the end of the simulation period, relative CMB deviations from simulation c0 increase disproportionately with decreasing length of the coupling interval and reach up to 26% depending on the scenario (Table 2). A systematic, trend in development across the four RCPs is, however, not observable.

Fig. 13. Deviations between the annual glacier-wide climatic mass balances of simulations c1, c2 and c3 and those of the uncoupled simulation c0 for RCP 2.6, 4.5, 6.0 and 8.5 (cf. Fig. 7).

In the warmest scenario (RCP 8.5), the ice cap thins from an initial maximum of 400 m to only 205 m under c3 or 245 m under c0. Large areas in the outlets become ice-free (Figs 5 and 11). For the other RCPs, changes in ice-cap geometry are less severe and depend less on coupling interval as climate change is reduced. For RCP 2.6, the final maximum ice thickness lies between 330 and 340 m depending on the coupling interval; the overall glacier outline remains close to present-day (see the 50 m lines in Fig. 5). For RCP 4.5, changes to the glacier outline, especially in the southwest, can be observed with a maximum thickness between 300 and 310 m in the centre of the ice cap. For RCP 6.0, the final maximum thickness is between 260 and 280 m, and changes in the cross sections and also glacier outline are clearly dependent on the coupling interval (Figs 5 and 11). Table 2 also indicates the volume of the ice cap at the end of the simulation, relative to its initial volume. This ratio also depends strongly on the coupling interval; for example, in the case of RCP 8.5, between 19% and 30% of the ice volume remains by the end of the century depending on the coupling scenario.

In summary, the effect of the topographic feedback is clearly visible from the cross sections (Fig. 11): neglecting the feedback leads to an overestimation of the ice thickness over the whole ice cap, in the central area as well as the outlets. When neglecting the feedback, the absolute elevation changes increase the more marked the climate change. Although our two shortest coupling intervals lead to fairly similar results, an even shorter interval may be needed to make reliable predictions.

Vestfonna mass loss corresponds to a SLR equivalent of between 0.3 and almost 1 mm by the end of the century in all our simulations (Fig. 9; Table 2). This mass loss increases strongly with increasing RCP strength. Volume changes in the coupled simulations (c1 to c3) are higher for shorter coupling intervals and lead to ∼18% additional volume loss in the case of c3 relative to c0. No clear trend from one RCP to the other is visible; the relative volume change of RCP 6.0 is, for example, slightly higher than than those of RCP 4.5 and RCP 8.5. Hence, accounting for the topographic feedback is equally important for all climate scenarios. Over the course of the simulations, volume loss increases faster than linearly with the length of the coupling interval.

4.3. Towards full coupling

Taken together, this means that no, or too long a coupling interval leads to a distinct overestimation of specific CMB across the entire ice cap. In conjunction with this, the final surface elevations and thus the mass loss in the central parts as well as via the outlets of the ice cap are considerably underestimated (Figs 5, 9 and 11). This is in agreement with findings for the Greenland ice sheet from, for instance, Goelzer and others (2013), who found a 14–31% higher sea-level contribution during the next 200 years for coupled model runs, and Edwards and others (2014b), who found an additional sea-level contribution due to the topographic feedback at the end of the century of ∼5%. Thus even our shortest coupling interval (10 years in simulation c3) may not be sufficient to adequately reproduce real-world CMB evolution and ice-mass loss. Most likely, the ice cap would experience even more negative CMB alongside an even stronger and faster thinning with full, i.e. annual, two-way coupling between mass-balance and ice flow model.

Table 2 and Figure 10 suggest that differences in CMB and ice-volume change between uncoupled and two-way coupled simulations as a function of coupling interval can be adequately expressed as a quadratic function. With this assumption it becomes possible to estimate the annual, glacier-wide CMB and the associated ice volume loss in terms of SLR equivalents for fully two-way coupled simulations (C4). c1 to c3 lead to relative deviations of up to 17–26% in CMB at the end of the simulation period compared with the uncoupled simulation c0, and we estimate that in a fully two-way coupled simulation (C4) the relative deviations increase to 20–30% (Table 2). Thus, CMB projections that do not use full two-way coupling overestimate the glacier-wide CMB increasingly as the simulation period increases. In these estimations the impact of the applied RCP is not homogeneous. No clear trend from one RCP to another is visible. Similar observations can be made for the ice volume loss (Table 2). Simulations c1 to c3 are up to 14–17% different from simulation c0; the estimated difference for C4 is 17–19%. There is no simple relation between strength of RCP forcing and the estimation of C4.

Coupling affects CMB and ice volume loss by different amounts. This is because the glacier-wide CMB values represent annually updated calculations while the ice volume loss is a cumulative quantity that is summed up over the simulation period. Moreover, the glacier-wide CMB is a one-dimensional quantity that is averaged over a stepwise changing (according to the different coupling intervals) glacier extent, while the ice volume loss is a 3-D quantity that is forced by both CMB and nonlinear ice dynamics.

To conclude, when comparing the estimated ice volume loss of an ideal simulation C4 to the calculated losses of simulations c0 to c3 the underestimation of the ice volume loss increases considerably with increasing length of the coupling interval (Table 2). While the shortest coupling interval (10 years in simulation c3) produces a 1–3% underestimation, this percentage increases to 6–7% in simulation c2, to 11–12% in c3 and to 17–19% in the completely uncoupled simulation c0.

4.4. Alternative: use of a lapse rate and its shortcomings

Using coupled simulations with different time steps it is possible to approximate the outcome of a fully two-way coupled modelling architecture. However, even completing simulations c0 to c3 for a small ice cap remains technically challenging. Establishing a coupling interface between the models that have been developed in different contexts, and hence feature codes of completely different structures and/or languages, is not trivial. Therefore, we present the lapse rate approach (simulation lr) as a potentially feasible alternative that allows models to be run totally independently.

Using our lapse rate approach produces good-quality results for the first half of the century, and even longer for the milder RCPs (Figs 9 and 11; Table 2). From the ice-cap cross section (Fig. 11) lr systematically overestimates the ice thickness expected for full coupling (approximated here by c3). For RCP 2.6, lr produces a result close to c3 over the whole simulation time. For RCP 4.5 the volume loss starts to deviate from c3 only at ∼70 years into the simulation. The situation is fairly similar for RCP 6.0, but deviations from c3 start slightly earlier. For RCP 8.5, the final ice thickness lies between those for c2 and c1; the total volume loss is only slightly higher than given by c1 and starts to differ from c3 from 55 years of simulation time onwards.

Because of the success of our simple lapse rate approach, we are fairly confident that using a (more sophisticated) lapse rate approach would be a good workaround to compensate for the topographic feedback when full coupling between climate (CMB) and ice-sheet model is not feasible. We now discuss what kind of improvements should be considered. The final quality of a lapse rate approach depends on its validity for each CMB component, the relative importance of this component for the mass budget, but also the sensitivity to topographic feedback of that component in the context of model coupling. Additionally, lapse rates may also vary over time or be different for different regions of the ice cap, such as in the ablation and accumulation areas.

Lapse rates as well as corresponding R ² values for the different CMB components are listed for the beginning and the end of the simulations in Table 3. These components are the air-temperature part of the ablation (first summand in Eqn (5)), a_T , its radiation part (second summand in Eqn (5)), a_R , the accumulation (Eqn (4)), c, and refreezing, r. Global R ² values decrease over the simulation for the stronger climate changes from 98% to 92%. a_T is the component which can clearly be best parameterized by a lapse rate (R ² of 97–99%). For a_R , at high altitudes, and to a lesser extent at lower altitudes, deviations from linear behaviour can be observed and become more pronounced during the simulations, in particular for more pronounced climate changes (R ² as low as 74%). Similar, but relatively smaller, effects can been seen in c. Accounting for one-third of the total glacier-wide mass gain, r can clearly not be parameterized by a simple lapse rate, even if this is not openly visible from the R ² value. As r is dependent on ablation, accumulation and rainfall, it does not show a linear behaviour with altitude and needs to be represented by different lapse rates below and above the equilibrium line. For the other components no distinct difference between ablation and accumulation area is observed, as was also noted by Helsen and others (2012) for the Greenland ice sheet. A clearly limiting factor is the time dependency of the individual lapse rates. This effect is most important for a_T , which changes by a factor of 3.7 between the start and the end of the simulation in the case of RCP 8.5, followed by a_R , which nearly doubles in the case of RCP 4.5 during the same period. A much less pronounced time dependency can be observed for c and r.

Table 3. Per cent variance accounted for (R ²) and lapse rates, λ (in mm ice a⁻¹ m⁻¹), at the beginning of the century (‘init’) and at 2010 for each RCP for each separate CMB component (i.e. radiation-driven part of the ablation, a_R , temperature-driven part of the ablation, a_T , surface accumulation, c, and refreezing rate, r)

In all scenarios applied here, the CMB is clearly dominated by ablation (Möller and Schneider, 2015). This amplifies the limitations introduced by the poor validity of a lapse rate approach for a_R and the strong time dependency of a_T . a_R accounts for up to just over half of the total glacier-wide ablation, but global radiation is little affected from the feedbacks in the coupling. Hence, surprisingly, the most limiting factor of our lapse rate approach turns out to be the air-temperature-dependent part of the ablation, a_T , despite the fact that it displays a very close to linear elevation dependency. On the one hand, air temperature shows a strong elevation dependency, hence this component is very sensitive to topographic feedback. In addition, ablation is the dominant contribution to CMB. On the other hand, there is a pronounced time variation of the elevation dependency, that needs to be captured. a_T , as written in Eqn (5), seems to follow the simple linear elevation dependency of the surface temperature; this is, however, misleading since ablation only takes place on those days of the year that are warm enough. Hence, when calculating the yearly ablation, the original elevation dependency is blurred by summing only over warm enough days, i.e. over more days in lower areas (and cooler climate scenarios) than in higher areas (and warmer climate changes). The lapse rate suggested by Eqn (5) represents only an upper limit for the lapse rate specific to a_T (4.5 mm ice eq. a⁻¹ m⁻¹), which would be reached if ablation took place over the whole ice cap every single day. Vizcaino and others (2010) also identify the changes in surface temperature and its elevation gradient as the main driver behind the topographic feedback between ice-sheet and mass-balance models.

In the absence of real model coupling, applying individual lapse rates with their time dependencies to the different components when correcting the uncoupled CMB fields could be considered as a more complex lapse rate approach. Furthermore, keeping a_R as in uncoupled model runs could be an option, since the assumption of a linear elevation dependency might introduce more error than neglecting the topographic feedback, which is small in the case of radiation. For r, replacing the linear dependency by a stepwise linear function would lead to additional improvement.

5. Conclusions

Running a full-Stokes ice-sheet model over a century with four different RCPs uncoupled from a CMB model leads to substantial underestimations of sea-level rise and overestimations of glacier-wide CMB. The mass-balance/elevation feedback between a time-evolving ice-cap surface and a climatic mass-balance model (or even further, a climate model) should therefore not be neglected. Our coupling experiments allow an estimate for the volume change in an ideal fully two-way coupled simulation, which is under-predicted by ∼18% in uncoupled model runs. The glacier-wide CMB is overestimated by 20–30% in uncoupled model runs. The remaining ice volume is overestimated by up to 50% (in the case of RCP 8.5), so predictions of the disappearance of ice caps are also severely affected when neglecting model coupling.

Using a lapse rate to account for changes in the topography when calculating the CMB, instead of a full model coupling, leads to considerable improvements compared to running both models uncoupled, at least for simulations of a few decades and modest climate-change scenarios. The most limiting factor in this approach is the time dependency of such a lapse rate. Using different approximations for the ice temperature and different temporarily fixed basal drag parameter fields also leads to variations in the predicted surface elevations and volume changes. However, the obtained results regarding the elevation feedback show no significant response to such variations. Changes in the basal sliding parameter, β, introduce locally important thickness variations although the total volume change is little affected by such changes. The more marked the climate changes, the less important are the effects of the different applied fields.

These results are qualitatively also valid for other ice caps and ice sheets, for which, because of their size, the impact of the underestimation of sea-level rise might represent a real concern. Hence, we strongly recommend that for reliable future sea-level rise predictions climate or mass-balance models and ice-flow models should be run coupled with using sufficiently short coupling intervals. Alternatively, full coupling can be replaced by deepened sensitivity studies, for example by a high-quality lapse rate approach or by elaborated extrapolations from longer to shorter coupling intervals.