2.1. Research location

The European Arctic archipelago of Svalbard has ∼34 560 km² of glaciated area (Reference Moholdt, Hagen, Eiken and SchulerMoholdt 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²; Reference Moholdt, Nuth, Hagen and KohlerMoholdt and others, 2010b) and Vestfonna (2340 km²; Reference BraunBraun 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 (Reference BraunBraun 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 (Reference Petterson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPetterson 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 (Reference BraunBraun and others, 2011; Reference Petterson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPohjola 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 (Reference BraunBraun and others, 2011; Reference Petterson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPohjola and others, 2011), suggesting a radical change in flow-velocity regime but with no clear indication of a typical surge (Reference SchäferSchäfer and others, 2014). However, past surges have been reported for Franklinbreen and other outlets of Vestfonna (Reference Błaszczyk, Jania and HagenBł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 (Reference MöllerMö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 (Reference Sauter, Möller, Finkelnburg, Grabiec, Scherer and SchneiderMöller and others, 2013). Accumulation shows high temporal variability, with changes of up to 100% between individual years (Reference BeaudonBeaudon and others, 2011). Spatial variability across the ice cap is governed by terrain elevation and snowdrift (Reference MöllerMöller and others, 2011b; Reference Sauter, Möller, Finkelnburg, Grabiec, Scherer and SchneiderSauter and others, 2013).

Numerical modelling studies by Möller and others (Reference Möller2011a, Reference Sauter, Möller, Finkelnburg, Grabiec, Scherer and Schneider2013) 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 (Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMö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 (Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMöller and others, 2013). This is consistent with remote-sensing-derived geodetic mass balance of the ice cap (Reference Nuth, Moholdt, Kohler, Hagen and KaabMoholdt and others, 2010b; Reference Nuth, Moholdt, Kohler, Hagen and KaabNuth and others, 2010).

2.2. Data basis and preparation

The study requires topographic data of VSF, i.e. bedrock and surface elevations. Bedrock elevation (Reference Petterson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPetterson 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 (Reference JakobssonJakobsson and others, 2008) for surrounding sea-floor. Details of the surface velocity data used for inferring the basal friction parameters are provided by Reference SchäferSchä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 (Reference MossMoss and others, 2010) are used (Fig. 2). These data are taken from the Earth System Model MIROC-ESM (Reference WatanabeWatanabe 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 Reference Möller and SchneiderMö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 Reference Möller and SchneiderMö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 Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMöller and others (2013). Air temperature downscaling was based on a combination of multiple linear regression and variance-inflation techniques (Reference Karl, Wang, Schlesinger, Knight and PortmanKarl and others, 1990; Reference HuthHuth, 1999; Reference Von Storchvon Storch, 1999; Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMöller and others, 2013). Precipitation downscaling combined multiple linear regression and local scaling techniques (Reference Widmann, Bretherton and SalathéWidmann and others, 2003; Reference Möller and SchneiderMö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 Reference MöllerMö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 (Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMö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.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 (Reference GagliardiniGagliardini and others, 2013). Here we use the methods described by Schäfer and others (Reference Schäfer2012, Reference Schäfer2014).

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 (Reference GlenGlen, 1952; Reference NyeNye, 1952), linking deviatoric stresses and strain rates. The temperature dependency of the viscosity is described by an Arrhenius law (Reference SchäferSchä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 (Reference Greve and BlatterGreve 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 Reference Arthern and GudmundssonArthern and Gudmundsson (2010) as detailed in Reference SchäferSchä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 (Reference SchäferSchä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 Reference SchäferSchäfer and others (2014) we use an anisotropic mesh utilizing the fully automatic, adaptive, isotropic surface re-meshing procedure YAMS (Reference FreyFrey, 2001) for the ice flow model. An initially homogeneously spaced two-dimensional footprint mesh was established using Gmsh (Reference Geuzaine and RemacleGeuzaine 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 (Reference Petterson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPohjola 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 (Reference CogleyCogley and others, 2011), is modelled using the temperature–net-radiation index model of Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMöller and others (2013). This model builds on that developed by Reference MöllerMöller and others (2011a) but incorporates monthly updated surface albedo-field calculations according to Reference MöllerMö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⁻¹ (Reference MöllerMöller and others, 2011a). The empirical parameters f _P1 (0.99 ×10⁻⁵ m⁻²) and f _P2 (0.79 × 10⁻² m⁻¹) were calibrated by Reference MöllerMöller and others (2011a) using in situ snow-pit data from Reference MöllerMöller and others (2011b). The transition from liquid to solid precipitation between +2 and 0°C is realized using a hyperbolic function (Reference Möller, Schneider and KilianMö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 Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMö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 Reference MöllerMö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 (Reference MöllerMö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 (Reference MöllerMöller and others, 2011a). Refreezing, r(x, y, z), is realized by applying the P _max approach of Reference ReehReeh (1991). Here we use P _max = 0.9 (Reference Möller, Finkelnburg, Braun, Scherer and SchneiderMö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.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 (Reference EdwardsEdwards 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.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 (Reference GoelzerGoelzer 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 Reference SchäferSchä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, Reference GoelzerGoelzer and others (2013), who found a 14–31% higher sea-level contribution during the next 200 years for coupled model runs, and Reference EdwardsEdwards 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 Reference Helsen, Van de Wal, Van den Broeke, Van de Berg and OerlemansHelsen 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 (Reference Möller and SchneiderMö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. Reference Vizcaino, Mikolajewicz, Jungclaus and SchurgersVizcaino 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.