Assessment of heat sources on the control of fast flow of Vestfonna ice cap, Svalbard

Understanding the response of fast flowing ice streams or outlet glaciers to changing climate is crucial in order to make reliable projections of sea level change over the coming decades. Motion of fast outlet glaciers occurs largely through basal motion governed by physical processes at the glacier bed, which are not yet fully understood. Various subglacial mechanisms have been suggested for fast flow but common to most of the suggested processes is the requirement of presence of liquid water, and thus temperate conditions. We use a combination of modelling, field, and remote observations in order to study links between different heat sources, the thermal regime and basal sliding in fast flowing areas on Vestfonna ice cap. A special emphasis lies on Franklinbreen, a fast flowing outlet glacier which has been observed to accelerate recently. We use the ice flow model Elmer/Ice including a Weertman type sliding law and a Robin inverse method to infer basal friction parameters from observed surface velocities. Firn heating, i.e. latent heat release through percolation of melt water, is included in our model; its parameterisation is calibrated with the temperature record of a deep borehole. We found that strain heating is negligible, whereas friction heating is identified as one possible trigger for the onset of fast flow. Firn heating is a significant heat source in the central thick and slow flowing area of the ice cap and the essential driver behind the ongoing fast flow in all outlets. Our findings suggest a possible scenario of the onset and maintenance of fast flow on the Vestfonna ice cap based on thermal processes and emphasise the role of latent heat released through refreezing of percolating melt water for fast flow. However, these processes cannot yet be captured in a temporally evolving sliding law. In order to simulate correctly fast flowing outlet glaciers, ice flow models not only need to account fully for all heat sources, but also need to incorporate a sliding law that is not solely based on the basal temperature, but also on hydrology and/or sediment physics.


Introduction
Recent studies suggest that an important contribution to sea level rise over the coming decades will be cryospheric mass loss in the form of discharge from fast flowing ice streams or outlet glaciers (Meier et al., 2007;Moon et al., 2012;Jacob et al., 2012;Tidewater Glacier Workshop Report, 2013).Therefore, understanding the response of fast flowing features to changing climate is crucial in order to make reliable projections (Moore et al., 2011;Rignot et al., 2011;Shepherd et al., 2012;Dunse et al., 2012).Observations dating back several decades show multiple modes of fast ice flow behaviour including permanently fast flowing outlet glaciers or ice streams, networks of ice streams that switch between fast and slow flow (Boulton and Jones, 1979), pulsing glaciers (Mayo, 1978), short-term velocity variations of fast tidewater glaciers (Meier and Post, 1987), and surging glaciers (Kamb, 1987;Dowdeswell and Collin, 1990;Dunse et al., 2012;Gladstone et al., 2014;Sund et al., 2009Sund et al., , 2011Sund et al., , 2014) ) that show occasional massive accelerations of a factor 10 or higher compared to their quiescent state.The underlying processes behind this range of behaviours are not yet fully understood and need to be addressed.
The flow in fast flowing glaciers exceeds speeds that could be sustained only by internal deformation -even at ice temperatures close to the pressure melting point, when internal deformation is highest.Fast glacier flow is therefore considered to be caused by basal motion through a combination of non-zero ice velocity at the bed, sliding over the bed and fast deformation of soft basal ice or subglacial sediments.Whereas ice deformation is relatively well explained today (Paterson, 1994), the physical processes controlling basal motion remain to be better understood.Many processes and feedbacks have been suggested to influence basal motion, including sub-glacial hydrology (Kamb, 1987;Vaughan et al., 2008;Bougamont et al., 2011;van der Wel et al., 2013), deformation of sub-glacial sediments (Truffer et al., 2000), heat production from sliding (i.e.friction heating, Fowler et al., 2001;Price et al., 2008), strain heating (Clarke et al., 1977;Pohjola and Hedfors, 2003;Schoof, 2004) or thermal instabilities (Murray et al., 2000).Common to most of the suggested processes is the idea that basal motion requires the presence of liquid water, and thus temperate conditions, at the base of the glacier.Hence, in order to understand the mechanisms of fast flowing ice, it is essential to study the processes maintaining and causing temperate basal conditions, as well as more generally the mechanisms leading to changes in the thermal conditions in the ice.For this purpose we examine different heat sources and their impact on basal temperatures together with the redistribution of heat through advection over the ice cap.
Long-timescale oscillatory behaviour of fast flowing ice streams can be solely explained by coupled flow and temperature evolution (Payne and Dongelmans, 1997;Hindmarsh, 2009;van Pelt and Oerlemans, 2012).However, shorter timescale oscillations, such as surges, require additional feedbacks or alternative mechanisms (Fowler et al., 2010).The mechanisms behind surges are poorly known, but have been suggested to be mainly hydrologically controlled on temperate (Alaskan type) glaciers (Kamb et al., 1985) or associated with thermal instabilities of polythermal (Svalbard type) glaciers (Clarke et al., 1984;Payne and Dongelmans, 1997;Murray et al., 2000;Fowler et al., 2001Fowler et al., , 2010)).The assessment of the heat sources that contribute to the basal thermal regime is also essential for the understanding of such temporal and spatial oscillations.
Here, we use a combination of modelling, field, and remote observations in order to study links between the thermal regime, heat sources and basal sliding in fast flowing areas on Vestfonna ice cap.Vestfonna is one of the two major ice caps of Nordaustlandet, the second largest island of Svalbard.Compared to the neighbouring Austfonna ice cap, Vestfonna has a relatively small area of thick slow moving interior ice.It is instead dominated by fast flowing outlet glaciers that extend from the coast to regions close to the ice divide.All of the fast outlet glaciers on Vestfonna are thought to be topographically controlled (Pohjola et al., 2011) and several of them are believed to have surge-type behaviour (Dowdeswell and Collin, 1990).We will especially focus on the recent speed-up of Franklinbreen (Pohjola et al., 2011), one of the major outlet glaciers on Vestfonna, and simulate ice flow with the Elmer/Ice finite element Full Stokes ice dynamic model (Gagliardini et al., 2013).In the context of ice dynamic models, the term basal sliding is used rather than basal motion, because basal motion is typically modelled by sliding of the ice over the bed.The in situ complexity is represented in the sliding law by a few free parameters (usually sliding or friction coefficients) -as long as specific sediment or hydrological models are not used.
While many large-scale flow models use spatially uniform parameters for the friction law, but only allow sliding when ice at the bed reaches the pressure melting point (Greve, 1997;Ritz et al., 2001;Quiquet et al., 2013), we solve an inverse problem to constrain spatially varying friction law parameters by determining the best match between model and observed surface velocities.This approach allows quantification of the basal sliding velocity, which can help to constrain the in situ processes (Morlighem et al., 2010;Pralong and Gudmundsson, 2011;Jay-Allemand et al., 2011;Habermann et al., 2012Habermann et al., , 2013)).It also improves the accuracy of the simulation-reproduced spatial patterns of observed surface velocities.It does not, however, give a direct relationship between temperature and in situ friction parameter.No predictions for the basal friction parameters at other times than the inversion are possible.From the model we extract information about basal frictional heating (due to sliding at the bed) and strain heating (due to internal ice deformation) as well as their possible evolution during the acceleration of Franklinbreen between 1995 and 2008.In addition, we use the Wright P max formulation (Wright et al., 2007) for refreezing to assess the role of latent heat release due to refreezing of percolating surface melt water in the snow pack (firn heating) that is advected through the ice.This allows us to identify the heat sources responsible for a temperate bed in the fast-flowing outlet glaciers.Comparing their time evolution provides insights in the driving mechanisms behind the observed recent acceleration of Franklinbreen and the conservation of fast flow in all outlets.
The research area and observational data have been presented in detail by Schäfer et al. (2012).Key features and additional data are described in Sect. 2 and Sect. 3 describes the ice flow model.In Sect. 4 we outline our different simulations.Results with respect to the dominant heat sources, relationships and feedbacks between fast flow, acceleration and the thermal regime are discussed in Sect. 5 before we conclude in Sect.6.

Research area and observational data
Vestfonna (VSF) is characterised by a varied surface topography with two main ridges and strongly pronounced fastoutlet glaciers (Dowdeswell and Collin, 1990, Fig. 1).In common with most glaciers in Svalbard, VSF is polythermal (Schytt, 1964;Palosuo, 1987).VSF was the target of a recent International Polar Year (IPY) project (Pohjola et al., 2011).The observational record extends back to the International Geophysical Year (IGY) 1957-1958(Schytt, 1964)), when data on surface elevation (using barometer methods) and ice thickness (from seismic surveys) were gathered.One focus of this work is Franklinbreen, the largest outlet glacier, which is located on the northwestern side of the ice cap and has recently accelerated (Pohjola et al., 2011;Braun et al., 2011).

Digital elevation models of surface and bedrock topography
For surface elevations 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 (Fig. 1).This DEM is completed with the International Bathymetric Chart of the Arctic Ocean (IBCAO, Jakobsson et al., 2008) for surrounding sea floor.A comparison between this NPI DEM (1990) and the SPOT-Spirit DEM (2007) shows a difference of less than 10 m on average over Franklinbreen (M.Braun, personal communication, 2010).This misfit is within the uncertainty of the two DEMs and hence the same DEM can be used for all our simulations spanning the period 1995-2011.The bedrock data are a combination of groundbased impulse radar and airborne radio-echo soundings (Pettersson et al., 2011, Fig. 1).The ground-based radar was deployed in the central part, while the airborne radar covers outlet glaciers and frontal areas.The interpolated combined DEM has a resolution of 500 m and a vertical accuracy of 25 m; 55 % of the total area of the interpolated DEM is within spatial autocorrelation of more than 0.6 and is deemed to be of satisfactory quality (for further details see Pettersson et al., 2011).Despite a partial inaccuracy due to the relative sparse bedrock data, the first-order variability in the bed topography, i.e. major subglacial valleys and peaks, are captured correctly, which is most relevant for modelling the whole ice cap.Exceptions are the southwestern tip (Idunbreen), as well as the lower parts of Bodleybreen and Rijpbreen, where errors in the DEM may be significant, as discussed further in Sects.4.1.2and 5.3.1.

Remote sensing data of surface velocities
The inversion modelling technique used to derive basal friction parameters, requires input of measured horizontal surface velocities.Tandem Phase ERS-1/2 1 day SAR scenes were acquired between December 1995 and January 1996.Surface ice velocities were calculated using SAR interferometry (InSAR) over most of VSF, apart from small areas over the lower part of the fast outlet glaciers, where dualazimuth offset-tracking was employed (henceforth "1995 velocities", Pohjola et al., 2011;Schäfer et al., 2012).Four ALOS PALSAR scenes were acquired between January 2008 and March 2008 with a 46-day time interval and velocities calculated using offset-tracking (henceforth "2008(henceforth " velocities", Pohjola et al., 2011)).For 2011, an ERS-2 SAR data stack, which was acquired in March/April with a 3-day time interval and processed with a combined InSAR and tracking approach similar to the 1995 data (Pohjola et al., 2011(Pohjola et al., ), is used (referred as "2011 velocities", unpublished data) velocities", unpublished data).In all cases, the vertical components of the velocities have been neglected during the calculation of horizontal velocities.
The displacement error in the InSAR data is 2 cm, which corresponds to a velocity error of 7 m yr −1 for Tandem ERS-1/2 SAR data (1-day time interval) and 2 m yr −1 for 3-day ERS-2 SAR data (Dowdeswell et al., 2008).By considering a matching error estimate of 1/10th of a pixel, the precision of offset-tracking is about 10 m yr −1 for the 2008 ALOS PALSAR data separated by a temporal interval of 46 days (Pohjola et al., 2011).In the 2011 ERS-2 data set (Fig. 1c), dual-azimuth offset-tracking was considered in the northern part of Vestfonna for the fast-flowing glaciers and here the matching error is estimated to be about 35 m yr −1 ; in the southern part SAR data of only one orbit is available.Consequently InSAR processing could not be used, and the error of range-azimuth offset tracking is very large, of the order of 130 m yr −1 .The 2008 and 2011 data sets do not cover the ice cap completely, and data gaps have been filled by inter-polation and smoothing, except for the southwestern corner (a region of slow flow), where the 1995 data were used to fill a larger data gap (neglecting possible variations in Gimlebreen).
The surface velocities are presented in Fig. 1 (before interpolation and patching).We observe in all three data sets the two very different flow regimes: slow ice flow over the central area of the ice cap and high velocities in the outlet glaciers.Between 1995 and 2008 a net speed-up of at least 100 % (doubling of speed) in the Franklinbreen outlet can be seen (Pohjola et al., 2011), which stabilised during 2008 until 2011.The southern branch continued to accelerate slightly, while the northern branch decelerated (Fig. 1).Franklinbreen is the outlet glacier showing the biggest changes since 1995.It reaches speeds comparable to other fast-flowing outlet glaciers in 2008-2011, which, nevertheless, are modest compared to other Svalbard surging glaciers (Hagen et al., 1993).

Thermal boundary conditions
Different thermal boundary conditions are required in the model, one being surface or air temperature.Svalbard's climate has a maritime character with cooler summers and warmer winters than is typical at such a high latitude (Möller et al., 2011).Mean monthly air temperatures on VSF do not exceed +3 • C, and winter monthly means fall between −10 and −15 • C with minimum values of the order of −25 to −40 • C (Möller et al., 2011).A lapse rate approach is used in the current study to prescribe the unaltered surface temperature (1) at the surface elevation S(x).We use a lapse rate γ = 0.004 K m −1 (Wadham and Nuttall, 2002;Wadham et al., 2006;Schuler et al., 2007).This value is close to the one adopted in other studies: γ = 0.0044 K m −1 (Pohjola et al., 2002).Liljequist (1993) found a slightly larger lapse rate of γ = 0.005 K m −1 from measurements between the summit of Vestfonna (known as Ahlmann summit) and the 1957/58 IGY station at Kinnvika.Data from the atmospheric model WRF (Skamarock et al., 2008;Hines et al., 2011) et al., 2007;Möller et al., 2011).We find a mean annual temperature of −7.7 • C and a mean winter temperature of −14.5 • C at sea level.
In addition, the geothermal heat flux is an important basal boundary condition.Contrary to Schäfer et al. (2012), who assumed a geothermal heat flux of 63 mW m −2 typical for post-Precambrian, non-orogenic tectonic regions Arrhenius law as in Paterson (1994), temperatures relative to pressure melting point (Lee, 1970), we take the value of 40 mW m −2 .This value is motivated by the measured gradients of profiles obtained by deep drilling on the Nordaustlandet ice caps (Zagorodnov et al., 1989;Ignatieva and Macheret, 1991;Motoyama et al., 2008).In the case of Nordaustlandet, ground surface temperature changes in the uppermost 1-2 km of the bedrock are most likely still influenced by the cold of the Weichselian period, explaining this lower measured value of 40 mW m −2 and leading to good simulations of an observed (via deep drilling) temperature profile on VSF (Motoyama et al., 2008), as explained in Sect.4.3.

Model description
The model equations are solved numerically with the Elmer/Ice model.It is based on the open-source multiphysics package Elmer developed at the CSC -IT Center for Science in Espoo, Finland, and uses the finite element method (Zwinger et al., 2007;Gagliardini and Zwinger, 2008;Gagliardini et al., 2013).More details on the model implementation for VSF can be found in Schäfer et al. (2012).Numerical parameters used in our study are summarised in Table 1.

Forward model
The ice is modelled as a non-linear viscous incompressible fluid flowing under gravity over a rigid bedrock.The force balance (quasi-static equilibrium) is expressed by the Stokes equations.The rheology of the ice is described by Glen's law assuming isotropic behaviour.The temperature dependency of the deformation rate factor, A(T , T pm ), is described by the Arrhenius law (Paterson, 1994) with parameters as in Table 1, where T is the temperature and T pm the pressure melting point.The evolution of free surface S is governed by a kinematic boundary condition using the climatic mass balance (Cogley et al., 2011) as input.S is assumed to be a stress-free surface, i.e. τ • n = 0, where n is the normal unit vector and τ the stress tensor.
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: where β is the basal friction parameter, t a unit vector in the tangential plane aligned with the basal shear stress, v || and τ || are the components of the velocity, v, and stress parallel to the bed at the base.We assume zero basal melting (v •n = 0).The basal friction parameter field β(x, y) will be inferred in this study from surface velocities using an inverse method (Sect.3.3).
On the lateral boundaries the normal stress is set to the water pressure exerted by the ocean, p w , for elevations below sea level, else we assume stress-free condition.Schäfer et al. (2012) use a temporally fixed depth-dependent temperature profile (here referred to as depth-dependent temperature profile)

Temperature model
where q geo = 40 mW m −2 is the geothermal heatflux, κ = 2.072 W K −1 m −1 the for the temperature range representative heat conductivity of ice and D(x) the ice depth (vertical distance to the surface at a given location x in the ice body).Unlike Eq. ( 3), in this work the temporal evolution of the temperature field is governed by the heat transfer equation, which reads where ρ is the ice density and Q a volumetric heat source.
The weak formulation of Eq. ( 4) constitutes a so-called variational inequality.A numerically consistent method to solve this inequality, which is implemented in our model, is explained in (Zwinger et al., 2007;Gagliardini et al., 2013) proved to show the best convergence properties as initial condition.The heat capacity, c, and heat conductivity, κ, are functions of temperature (Ritz, 1987), turning Eq. ( 4) into a non-linear problem: At the upper boundary a Dirichlet condition is imposed on the temperature, T , using the parameterisation described by Eq. (1).At the bed a jump in the normal component of the imposed geothermal heat flux, q geo = κ gradT •n, is given by the surface production due to friction heat, q f = v || τ || , which in the case of a temperate base is completed by the latent heat sink caused by the basal melt rate.
The volumetric heat source, Q, comprises strain heat, Q s = 2µ , where is the second invariant of the strain rate tensor and µ the effective ice viscosity, as well as latent heat, Q l , released during refreezing of percolating melt water inside the firn layer (firn heating).Q l is calculated using the P max model of Wright et al. (2007) as used by Zwinger and Moore (2009).P max is defined as the maximum proportion of the annual snowfall that can be retained by refreezing before runoff occurs and was originally chosen to be 0.6 (Reijmer et al., 2012).The depth-integrated amount of energy Q l , P max and the annual snowfall, B, are related through the latent heat of fusion, L, by The model of Wright et al. (2007) does not use such a constant P max value, but provides a simple, yet realistic method of calculating P max as a function of the mean annual and mean winter temperature and provides simultaneously an estimate for the firn heating Q l .The energy to warm the uppermost part of a glacier from the end of the winter to postrefreezing temperatures is estimated and identified with the heat sink available to be filled with latent heat, Q l .Following this approach, different characteristic shapes of the timeaveraged temperature-depth profiles, T (d), in summer and winter are used (Wright et al., 2007) to calculate the latent heat at each point and each depth of the glacier: where d is the depth below the surface.There are three free parameters in this firn heating parameterisation: T a and T w are the annual and winter mean air temperatures respectively set according to Sect.2.3; d ice is the typical penetration depth of the annual temperature cycle, which is kept as a free parameter and tuned to reproduce the measurements in the deep ice core (Motoyama et al., 2008), see Sect.4.3.
The resulting non-uniform volume heat source is deduced by the difference of internal energy defined by the difference T (d) between the seasonal profiles Eqs. ( 8) and ( 9).It decreases steadily from the surface to the penetration depth.In all layers below d ice it is zero: P max can then be determined by Eq. ( 7).

Inverse model
A variational inverse method (Arthern and Gudmundsson, 2010) is used in this study to infer the spatially varying basal friction parameter β(x, y).This approach is similar to Schäfer et al. ( 2012), but with the addition of a regularisation term (Morlighem et al., 2010;Habermann et al., 2012).
The method iteratively applies a Neumann and a Dirichlet condition at the upper free surface as introduced by Maxwell et al. (2008).In the Dirichlet problem the Neumann free upper surface condition is replaced by a Dirichlet condition, where the observed surface horizontal velocities are imposed: where v hor (x) and v obs (x) stand respectively for the modelled and observed horizontal surface velocities.In zdirection, (τ • n) • e z = 0 is imposed on S, where e z is the unit vector along the vertical.To avoid unphysical negative values, the friction parameter field, β(x, y), is expressed as β = 10 α and the minimisation of the cost function is performed with respect to its logarithm, α.The cost function, which expresses the mismatch between the two solutions for the velocity field with different boundary conditions on the upper surface, S, is given by where the superscripts "N" and "D" refer to the solutions of the Neumann and Dirichlet problems, respectively.To avoid unphysical small wavelength variations in α and to ensure to find a stable unique solution, a Tikhonov regularisation term, J reg , penalising the first spatial derivatives of α is added to the total cost function, J tot : The Cryosphere, 8, 1951-1973, 2014 www.the-cryosphere.net/8/1951/2014/where λ is a positive parameter (see Sect. 4.1.1 for its choice).The minimisation of the cost function is thus a compromise between best fit to observations and smoothness of α, determined by the tuning of λ.The minimisation algorithm is described by Gillet-Chaulet et al. (2012) and Gagliardini et al. (2013).

Meshing
Anisotropic mesh refinement is now increasingly used in numerical modelling, especially with the finite element method since the mesh resolution is a critical factor.Schäfer et al. (2012) have investigated effects of varying the resolution in the context of this inverse method.Here we use again the mesh established with the fully automatic, adaptive, isotropic surface-remeshing procedure Yams (Frey, 2001).A 2-D footprint mesh was established according to the glacier outline on the 1990 NPI map and adapted using the metric based on the Hessian matrix of the observed 1995 surface velocities.Horizontal resolution varies between 250 and 2500 m.Finally the mesh was extruded vertically in 10 equidistant terrainfollowing layers according to the bedrock and surface data.
The obtained mesh consists of linear wedge type and hexahedral prism elements.In the simulations involving firn heating, the mesh was extruded vertically in 20 layers.The lower 10 layers have the same thickness, while in the upper 10 layers thickness is reduced towards the surface layer following a power law.The robustness of the total vertical layer number has already been verified (Schäfer et al., 2012), and doubling the number of boundary layers also leads to robust results in the runs including firn heating.

Simulations
In this section we present the set-up of our simulations.Four types of simulations are conducted and summarised in Table 2: -Simulations dealing with the thermo-mechanical spinup.These serve as starting points for all other simulations (Sect.4.1).Since an ideal thermo-mechanical spin-up is unfeasible, we describe our alternative initialisation.First, the distribution of the basal friction parameter regulating the velocity field is determined.Then, a purely mechanical spin-up is conducted followed by the calculation of a temperature field making certain assumptions.
-Simulations to investigate the importance and influence of the mechanical heating (strain and friction heating) (Sect.4.2).
-Simulations to calibrate our firn heating formulation (Sect.4.3) with a measured temperature record in a deep borehole.
-Short prognostic simulations (Sect.4.4) aiming to reproduce the observed acceleration of Franklinbreen and identify the underlying driving mechanisms.

System initialisation
Starting from a DEM purely based on observed data raises an inherent problem since a consistent initial condition for the thermo-mechanically coupled system is required.Thermal initial conditions are critical in modelling of polythermal glaciers or ice sheets because of the energy storage capacity of ice, the low advection/diffusion rates on the glacier and the strong thermo-mechanical coupling via the ice viscosity.
An ideal spin-up would demand a transient run starting from deglaciated conditions with a long enough spin-up time requiring realistic forcing (temperature, mass balance), as well as knowledge of the strongly time-varying velocity field.Air temperature and precipitation records might exist over a long enough time, however the temperature distribution at a given instant is driven by the past evolution of advection, diffusion of heat and heat sources, and hence by the past velocity field.
In the absence of such ideal spin-up, we decouple inversion for basal friction parameter from the thermo-mechanical problem as detailed in the following sections.

Inverse simulations to derive spatial patterns of the basal friction
The method of Schäfer et al. ( 2012) is followed with some improvements: correct marine boundary conditions are applied, a regularisation term in the cost function is added and a better minimisation algorithm is used.The best value of the regularisation parameter, λ, in Eq. ( 13) is determined by L curve analysis (Hansen, 2001) from a plot displaying J reg (smoothness of the friction parameter) as a function of J 0 (match to observations).This analysis is done for simplicity only once using the 1995 velocity data set.For this purpose, we use a fixed temperature distribution given by the depth-dependent profile (Sect.3.2), as done by Schäfer et al. (2012).This is justified since the result of the inversion depends only very little on the small thermally induced variations in the ice viscosity, as shown in Schäfer et al. (2012).We find that J 0 is minimised by setting λ = 10 5.0 (unit system MPa m yr), which also leads to acceptable smoothness in β (simulation beta95).Modelled and observed velocities show a close match, as in Schäfer et al. (2012).
In a similar way, we conduct inversions for the basal friction parameter with the 2008 and 2011 surface velocity fields (Fig. 2, simulations beta08 and beta11).These runs are conducted also with the 1990 surface DEM, since no complete additional surface DEM is available.Nuth et al. (2010) and Moholdt et al. (2010) 13), log scale, scaled by 10 9 -when iterating inverse method and temperature steady-state calculation (1995 velocities, simulation beta95T).The graph is cut at 1.6 × 10 9 for better visibility, higher values in the first iterations are thus not visible.In orange (void circles), the evolution of the cost function when not coupling with temperature is shown.
VSF, with local values up to 1 m yr −1 in the south.It has been shown (Schäfer et al., 2012) that surface variations of this order (or higher) between 1995 and 2008 or 2011 do not significantly affect the friction parameter fields derived from the inverse method.Some error is, however, expected from changes in the surface slope resulting from these small surface elevation changes (Joughin et al., 2004).The same value for the λ parameter as well as the same inhomogeneous mesh have been used, the latter to facilitate comparison.
An iteration scheme between inversion and temperature calculation has been tested (simulations beta95T and beta08T).The depth-dependent temperature profile (Eq. 3) was used in the first inversion for β, then steady-state calculations of the temperature field (accounting for friction and strain heating) and inversions were run alternately.The resulting β distribution reveals small changes compared to keeping the temperature fixed to the depth-dependent profile, showing a certain robustness of β towards changes in temperature.Nevertheless, the value of the cost function has been decreased with the iterative scheme (Fig. 3), in line with improved match between observed and computed surface velocities.Convergence of this iteration was assessed through the cost function and stopped once the cost function stabilised.Convergence of the steady-state temperature field was ascertained through visual inspection.
The effect of firn heating on the resulting β distribution has been studied separately (simulation beta95Tfirn).This is motivated both by the need to save computing resource (firn heating simulations require higher vertical resolution), and because of the greater uncertainties associated with this heat source compared to the others.It causes very little change to β.For all further simulations in the current study, the distribution of β obtained with the iteration scheme (but with firn heating omitted) is used (Fig. 2a, iteration No. 120 in Fig. 3).

Surface relaxation
Remaining uncertainties in the model initial conditions (including uncertainties in the model parameters, as well as the domain geometry), lead to ice flux divergence anomalies (Zwinger and Moore, 2009;Seroussi et al., 2011), resulting in a non-smooth vertical velocity field.Because of the importance of the vertical velocity field for advection of cold ice from the surface, and to smooth out these ice flux divergence anomalies, the free surface is relaxed before being used for further simulations.The relaxation simulations lasted for three years under mean present-day climatic surface mass-balance (Möller et al., 2011), simulations surfre-lax95 and surfrelax08.They were initialised with output from the inversion-temperature iterations beta95T, beta08T.
A short time step (0.1 year) was chosen to guarantee temporal resolution of artificially strong surface changes induced by the remaining uncertainties (Zwinger and Moore, 2009).
Visual inspection of the smoothness and magnitude (Fig. 4) of the vertical velocity field as well as surface elevation adjustments were used to determine the end of the relaxation procedure.The largest changes to the mesh (Fig. 5) occur in the southwestern corner and in some outlet glaciers, where there is a significant paucity of bedrock radar data (see Pettersson et al., 2011, for radar coverage).Some other less important changes are visible in northeastern VSF -again in areas with sparsely covered bedrock data.
A more complex spin-up scheme involving an iteration between surface relaxation, inverse method and temperature calculation was also tested for a single combination of surface velocity data and included heat sources in the temperature calculation (simulation surfrelax95c).This procedure  requires huge computational effort and does not lead to visible improvements in the results (β field, temperature field, and surface corrections) and is hence not used in this work for several different combinations of surface velocity data (three possibilities) and included heat sources (six possibilities).

Thermal initialisation
In the absence of an ideal thermo-mechanical spin-up, a steady-state temperature field was computed after the mechanical relaxation, even though non-steady-state conditions have occurred on VSF between 1995 and 2008, and probably on longer timescales.We found characteristic timescales to reach such a steady state to be of the order of several hundreds of years (not shown in this paper).This will lead to over-or underestimations of the temperature depending on the past state of each outlet glacier.Similarly, Seroussi et al. (2013) addressed the question of thermal initial conditions on the Greenland ice sheet and came to the conclusion that steady-state temperatures based on present-day conditions are, nevertheless, a reasonably good approximation, both for calculations of basal conditions and century-scale transient simulations.

Temperature steady states including mechanical heating
Strain and friction heating (Sect.3.2) are effective locally in the outlets and given by the mechanical model.To discuss their influence on the thermal regime of the ice cap, temperature steady states for various combinations of the heat sources are calculated (Fig. 6, simulations 1995ssA, 1995ssSH, 1995ssFH, 1995ss, 2008ss, Table 2).In these simulations only the internal ice temperature is allowed to evolve; the surface temperature, velocity field and geome-try are kept fixed.Simulated temperature-depth profiles are extracted (Fig. 7) in chosen locations that are indicated in Fig. 6.These are two locations of well-surveyed measurements sites on Franklinbreen and Frazerbreen, in the ablation and in accumulation areas respectively; as well as the location of the drill hole.

Calibration of the firn heating formulation
In order to calibrate the parameters (penetration depth, mean annual and winter temperatures) in the applied firn heating parameterisation (see Sect. 3.2), we study the temperature profile at the location of the drill hole (see Fig. 6 for the exact location).Using only strain and friction heating (simulation 1995ss), the measured temperature profile (Motoyama et al., 2008(Motoyama et al., , recorded in 1995) ) cannot be reproduced, even close to the bedrock, where these heat sources are most effective (Fig. 8, green line compared to the data in red).
In our firn heating formulation we assume that the penetration depth increases linearly from 0 m at the elevation of an average firn line to the maximum penetration depth d ice at the summit, leading to the effect that firn heating increases with altitude above the firn line.This is also in line with the usual approach of calculating refreezing as a fraction of winter accumulation, which is best described on VSF by an elevation gradient (Möller et al., 2011).In reality the melting should be largest at low (warmer) elevation, but this is at least partially counterbalanced by the formation of ice lenses inhibiting penetration of melt water.Ice lenses are defined as discrete anomalies in density of the firn column observed by density measurements and by geophysical scanning (DEP) and by ocular inspection of the ice facies.The latter separate the ice facies due to difference in void space/air bubble content.These observations are standard when analysing ice cores.
The Cryosphere, 8, 1951-1973, 2014   2) profile including advection, but no additional heat sources (green) (1995ssA), (3) profile including advection and strain heating (blue) (1995ssSH), (4) profile including advection and friction heating (pink) (1995ssFH), (5) profile including advection and both strain and friction heating (light blue) (1995ss), ( 6) profile including also firn heating (orange), which is identical for the equilibrium state firnss and the following time-evolving simulation firnevol.The locations of these profiles are indicated in Fig. 6.Ice lenses are more effective with more melt, thus in lower altitude.Assuming increasing firn heating with increasing firn thickness (altitude) is hence a simplification of these competing effects.The mean elevation of the firn line was digitised from several satellite pictures: Landsat July 1976, September 1988 andAugust 2006;Spot July 1991 andAugust 2008;Aster August 2000 andJuly 2005.Two of these lines (August 2008 and September 1988) have been excluded since the firn lines are located at exceptionally low elevations probably due to abnormally early fresh snow.We observe little change over recent decades, as found by Möller et al. (2013).Since the firn line elevation is approximately uniform over most of the ice cap, a single mean elevation for the firn line is assumed ignoring any other spatial variations.We estimate this mean elevation to be 410 m a.s.l.(Fig. 9).This is consistent with estimates of the average equilibrium line elevation (326 m a.s.l., Möller et al., 2013), since on Svalbard glaciers the equilibrium line is typically located significantly lower due to extensive superimposed ice formation (Möller et al., 2011;van Pelt et al., 2012).For future prognostics with climatically varying mass balance the firn line elevation could be parameterised by the elevation of the equilibrium line.
Using the mechanical set-up as described in Sect.4.2, a variety of simulations with different parameters of the firn of the firn heating is added or not.The position of the ice core is indicated in red, the position of the two locations used in Fig. 7 in white.The firn line is drawn in purple, the 20 m ice thickness isoline in light grey.The latter is a good indicator for underestimated ice thickness in areas where fast sliding seems to be in strong contradiction with too cold basal temperatures.
heating formulation were conducted.Different initial conditions, steady-state and time-evolving simulations have been tested and lead to the following conclusions: 1.The measured inflection in the upper part of the temperature profile is a transient effect occurring on decadal timescales.Fig. 8 illustrates the smoothing of the inflexion and its propagation towards the bedrock when approaching equilibrium.
2. Significant changes of temperature at the bedrock induced by firn heating occur in simulations over timescales long enough to transport heat to the bottom (centuries).Such simulations lead to temperature profiles similar to steady-state profiles (Fig. 8), but with a spatially uniform warm or cold shift.
Consequently we hypothesise that the measured borehole profile can be modelled by a succession of a steady-state followed by a time evolving simulation with different surface boundary conditions: low firn heating for the steadystate simulation and an increase in firn heating prescribed for the more recent time-evolving simulation, which has not yet reached equilibrium at present day.
In order to tune the model to match drill hole observations (Motoyama et al., 2008), we make the hypothesis that a boundary condition change can be represented in our firn heating parameterisation by the penetration depth parameter.We keep the surface temperatures fixed to the observations of the weather stations as stated earlier.A first run of the model (simulation firnss) to equilibrium temperature using a maximum penetration depth of 13.2 m leads to a good match with the measurements in the lower part of the drill hole, see Fig. 8 (blue line).A prognostic run (simulation firnevol) over 35 yr with an increased penetration depth of 17.6 m starting from this equilibrium state allows for a reasonably good match with the observed peak in the upper layers (black line).Because of advection, on a long enough timescale, firn heating affects the thermal regime of the whole ice cap including below the firn line.The horizontal distribution of the modelled temperature at the bedrock including firn heating is shown in Fig. 9 and the vertical in Fig. 7 as well as Fig. 8.

Prognostic simulations over the period 1995-2008
In addition to the steady-state experiments with the different heat sources, we conduct prognostic simulations with three different prescribed temporal evolutions of β to simulate the recent acceleration of Franklinbreen.We analyse the connected evolution of all system variables to get a better understanding of the underlying mechanisms.The three simulations are run with full thermo-mechanical coupling starting from the relaxed 1990 DEM (surfrelax95), forced by mean present-day climatic mass balance (Möller et al., 2011), and with constant surface temperature and glacier extent.The three basal friction parameter evolution scenarios are: 1.The basal friction parameter kept constant at the 1995 pattern (simulation const13).
2. A sudden switch after five years to the 2008 pattern (which differs from the 1995 pattern mainly by the acceleration of Franklinbreen, simulation sudden13).

A linear change from the 1995 to 2008 values (simulations lin13, lin13firnss and lin13firnevol).
All simulations span 1995-2008.Each simulation was run three times -excluding and including firn heating (firnss or firnevol) to study different scenarios.In the simulation excluding firn heating, temperature is initialised to the 1995 steady-state temperature profile 1995ss.In the simulations including firn heating, temperature is initialised to firnss and firnevol respectively.Changes in surface elevation and basal temperatures are shown in Figs. 10 and 11.

Implications of inferred basal friction parameter distributions
The basal friction parameter, β, is a crucial parameter in simulating the thermodynamical regime of VSF.It is a key control on sliding velocities, which govern both friction heating and heat advection.As already shown by Schäfer et al. (2012), the results of an inverse method used to derive the spatially varying basal friction parameter are largely unaffected by the temperature distribution (Fig. 2).This is because temperature, which impacts on deformation, does not feature in the sliding relation, and in VSF outlet glaciers sliding dominates over deformational velocities.Conversely, the temperature distribution shows a high sensitivity to the derived basal friction parameter, since even small changes can introduce important changes in the velocity field.Such changes are especially important for the vertical component, which then affects the heat redistribution through advection.Surface relaxation (Sect.4.1.2,simulations surfre-lax95,08) reduces this sensitivity by producing smoother velocity fields.
Variations in the basal friction parameter distribution across the three periods (Fig. 2) are due to large variations in observed surface velocities and indicate the importance of a time-evolving basal friction parameter based on the underlying physical processes.When comparing the obtained basal patterns from 1995, 2008 and 2011 (simulations beta95, beta08, beta11), the fine structure of the basal friction parameter in some of the outlet glaciers differs slightly.The most In (a) the texture is the 1995ss temperature distribution (temperatures relative to pressure melting).The light blue corresponds to the extent of the "sliding area" for 1995ss, dark blue for 2008ss and pink for the end of lin13.In (b) the temperature at the end of the 13 yr simulation when neglecting firn heating (lin13) (temperatures relative to pressure melting) is used as texture -same colour scale as (a).The light colours present the "sliding area" when neglecting firn heating (lin13), dark colours when including it (lin13firnss); blue lines represent the start of the 13 yr simulation and pink the end.
striking change remains the acceleration of Franklinbreen from 1995 to 2008 featured by a strong increase in basal sliding.The 2011 β-distribution mainly reflects the different changes in velocity pattern in the two branches of Franklinbreen: the northern one is decelerating, while the southern one continues to accelerate.In all outlet glaciers a distinct spatial variation of β can be seen, indicating that a sliding law also needs to reproduce these variations; especially since Schäfer et al. (2012) have shown that spatially constant friction parameters specific to each outlet glacier do not allow reproduction of the observed velocity structure within the outlet glaciers.A sliding law solely based on the presence of temperate ice at the base could not reproduce such a fine structure.

Interpretation of a temperature profile from deep drilling
As shown in our simulations concerning the calibration of the firn heating formulation (Sect.4.3, simulations firnss, firnevol), the observed shape of the temperature profile measured in 1995 in the borehole (Motoyama et al., 2008, Fig. 8) cannot be explained by an equilibrium temperature profile.Our model-supported interpretation requires a recent perturbation away from an earlier (close to) equilibrium state, caused by a change in the surface conditions.This can be motivated by the fact that ice cores elsewhere in Svalbard indicate that periods of firn heating and percolation have been frequent in the last 500-1000 yr (van de Wal et al., 2002;Divine et al., 2011).The proportion of ice lenses which indicate periods of near-zero ice surface temperatures increased from 33 % during the Little Ice Age to 55 % in the 20th century (Pohjola et al., 2002).
van de Wal et al. ( 2002) came to a similar conclusion when reconstructing the temperature record in the Lomonosovfonna Plateau (northeast of Billefjorden/Isfjorden, Spitzbergen).However, they kept the surface temperature as tuning parameter.Their obtained surface temperature is too high and induces a shift of a few Kelvin.Hence model and data fit well in the lower part, but the surface values are unrealistically warm.They conclude a change in surface conditions in the 1920s from their model and find confirmation for this by comparing to the mean air temperature record at Svalbard airport starting in 1910.Discrepancy between our model-implied change in the 1950s or 1960s and the actual climatic record can be explained by various facts: first, as stated earlier, the uncertainty in the basal friction parameter strongly impacts the evolution of the temperature distribution through advection.Second, only one data set of deep borehole temperatures is available for model calibration.Lastly, our approach might be too simplistic, especially with respect to the assumptions of spatial or elevation dependencies, and the use of penetration depth as the only calibration parameter.Surface temperature variations also lead to temperature variations at depth but with different timescales than the penetration depth and temperatures certainly have not been constant in the past.It is known for example that there were periods warmer than the first decade of the 21st century (D'Andrea et al., 2012).
The calibrated penetration depths (13.2 and 17.6 m, Sect.4.3) exceed the expected values, since measured densities (Motoyama et al., 2008) reach values over 850 kg m −3 at 10 m depth and below, i.e. values of ice or ice lenses with very slow percolation.These unrealistically high values for the penetration depth can be explained by the omission of firn layer compressibility in the mechanical model  (Zwinger et al., 2007).This implies that our approach should be considered as qualitative rather than quantitative.
With respect to this more recent change in the conditions on the surface (simulation firnevol), our model predicts that combined advective and diffusive processes will take over 100 yr to propagate this signal to the base in the centre of the ice cap (see Fig. 8) and over the outlets (not shown).Thus for studying basal processes, even for prognostic simulations up to a century, we can neglect the effects of this change in firn heating.Conversion of the latent energy released at the location of the ice core and comparing to snow fall corresponds at the location of the drill hole to a P max value of 0.9, which is in the expected range (Wright, 2005), increasing confidence in our model.
A discrepancy between our equilibrium profile and the data is also apparent in the middle of the depth profile.We give several possible explications for this: either the ice cap had not yet reached the first thermal equilibrium corresponding to the first penetration depth of 13.2 m before the change in surface boundary conditions occurred (see the 50/100/150 years etc. graphs in Fig. 8 for the shape of such profiles in the lower and middle part).Or other reasons could be our simplified assumption of a constant surface temperature, the impacts of uncertainties in advection (Sect.5.1) or the geothermal heat flux.With the model formulation representing latent heat release due to refreezing we are qualitatively able to reproduce the observed profile, indicating that we identified the driving mechanisms behind the measured distribution.

The role of heat sources for VSF fast-flowing outlet glaciers
A comparison between surface and sliding velocities at the bedrock (see Fig. 12a) clearly shows that sliding dominates the ice dynamics at the fast-flow areas of VSF, which even holds for increased deformation due to temperature-lowered viscosity (Schäfer et al., 2012).Schäfer et al. (2012) further showed that the temperature distribution has little impact on both surface velocities and basal friction parameter obtained with the inverse method.Therefore we focus on the impact of temperature and the respective heat sources on the onset and maintenance of fast flow.

Interpretation of calculated temperature steady states
Even though the ice cap had probably neither in 1995 nor in 2008 reached a temperature steady state, examination of the steady-state distributions gives insights in the possible impacts of the heat sources.

Strain and friction heating
Strain heat integrated over the whole ice column is of the same order of magnitude as friction heat at the bed (Fig. 12b).
It is mainly confined to the shear margins and corresponds well to the few areas of non-negligible deformational velocities (Fig. 12a).Friction heat is greatest at the transition of fast and slow flow, i.e. at the margins of the southern fastflowing outlets as well as in the areas of changing basal friction parameter of Franklinbreen (Fig. 12c).Simulation 2008ss shows that it is also important in areas of very high velocities at Franklinbreen (Fig. 13c).Larour et al. (2012) observed a similar spatial distribution of strain and friction heat.
In contrast to Pohjola and Hedfors (2003), who investigated fast flow in Antarctica using a one-dimensional numerical thermodynamic model, the impact of strain heating on the temperature field is found to be very small on VSF, see Figs. 6b and7, simulation 1995ssSH. Brinkerhoff et al. (2011) found for some Greenland outlets even less impact of strain heating.Friction heating allows temperatures close to pressure melting to be reached in a small area of the base of Franklinbreen (simulation 1995ssFH, Fig. 6c).In 2008ss (Fig. 6d) the temperature distribution shows colder areas, probably resulting from transport of colder ice from central regions towards the outlet.None of the southern outlets reaches a temperate base.
While friction heat might have played a role for the onset of the acceleration of Franklinbreen through a combined thickness-friction heat feedback, it cannot sustain fast flow in any of the outlets: in 1990 Franklinbreen was the thickest of the outlet glaciers and the only outlet glacier featuring some areas with sub-melt sliding (Echelmeyer et al., 1987;Hindmarsh and Le Meur, 2001), even in the absence of additional heat sources.By its reduced ice velocities it might have thickened enough to allow basal ice to approach pressure melting through insulation and in turn triggering sliding.Thickening of Franklinbreen especially in the lower part between 1990 and 2005 is confirmed by Nuth et al. (2010).However, there is no evidence in the calculated temperature distributions for friction heat being responsible for maintaining fast flow in any of the outlets.We also rule out a simple thickness feedback for the recent reduction in acceleration: no complete surface DEMs from different periods are available, though neither Nuth et al. (2010) nor Moholdt et al. (2010) observed a thinning on Franklinbreen between 1990-2005or between 2003-2008. Nuth et al. (2010) observe a balanced or slightly positive volume change over Franklinbreen (average of 0.06 ± 0.12 km 3 yr −1 ) for the period 1990-2005.Even though their result is subject to a large error, it seems unlikely that the observed recent reduction in accel-eration between 2008 and 2011 is driven by a mechanism involving thinning.

Firn heating due to melt water refreezing
Firn heating is important for the general thermal regime of the ice cap and has larger impacts than friction or strain heating in most of the regions (Figs. 7 and 9, simulation firnss).It is not only efficient above the firn line, but with a certain delay by advection also in all other parts of the ice cap.Advection timescales from the centre of the ice cap (region above the firn line) to the lower parts of the outlets are estimated to be of the order of centuries in our model (Fig. 14).Our model shows -assuming steady state -an increase in basal temperature in the onset areas as well as over the whole area of all fast flowing outlets (Fig. 9).There we observe a good correlation between temperate base and fast flow.
In the lower parts of some of the southern-facing outlet glaciers our calculated temperatures remain cold and we observe areas which are cold-based yet fast flowing (Fig. 9).This apparent contradiction of observed sliding over very cold bed is a consequence of our approach: basal friction parameters are inferred from observed surface velocities without direct coupling to temperature.Brinkerhoff et al. (2011) discussed the possibility of sliding over a cold base or underestimation of ice deformation due to neglecting ice anisotropy.Here, the apparent cold-based sliding occurs mainly in areas where the bedrock elevation is poorly  (Ahrens et al., 2005) assuming the basal friction parameter beta95T.The trajectories are coloured accordingly to the advection time to its lower end (blue colour).The firn line is depicted in orange, as texture the temperature field 1995ss is drawn.known, as confirmed by the 20 m ice thickness iso-line in the two problematic outlets (Fig. 9).Irregularities in the bedrock data and analysis of the 1995-2008 prognostic simulations (Sect.5.3.2) confirm this likely local underestimation of the ice thickness.Improvement of the bedrock data by some control method should be considered, as for example done by Morlighem et al. (2013) or van Pelt et al. (2013).

Simulating the acceleration period 1995-2008
By examining the temperature steady-state distributions, information about temporal evolution is disregarded.While the lower part of the temperature measured at the location of the deep drilling was in 1995 very close to the steady-state firnss, steady state has not necessarily been reached over the whole ice cap.We conduct simulations of the acceleration period to get further insights into the temporal evolution.
We first highlight some general observations of such prognostic simulations with different friction parameter evolution scenarios (Sect.4.4, Fig. 2) from the simulations without firn heating (const13, sudden13, lin13) focusing on Franklinbreen.Simulations sudden13 and lin13 exhibit a clear thinning of the onset area of Franklinbreen relative to const13 (Fig. 10).This is more pronounced (up to 25 m compared to 20 m) in simulation sudden13 resulting from the timeintegrated ice flux, which is higher the sooner the velocities increase.The pronounced increase in thickness at the terminus of Franklinbreen can be interpreted at least partly as a model feature caused by fixing the lateral extent of the ice cap and neglecting enhanced mass loss due to calving.Other outlet glaciers, especially Rijpbreen and Bodlebreen, are highly influenced by errors in the bedrock DEM and thus not discussed.
As expected, velocities, friction and strain heat remain constant during simulation const13.In simulation sudden13 a sudden jump after the change in the basal friction parameter is observed, while in simulation lin13 the variables change synchronously with changes in the basal friction parameter.The 2008 variables are very similar in simulations sudden13 and lin13.The velocities are not only faster than in 1995, but show the two pronounced branches of Franklinbreen.Friction heat values are fairly high up to the centre of the ice cap, also emphasising the developing two branches (Fig. 12 compared to Fig. 13).Strain heat increases at the lateral margins of Franklinbreen over the whole ice depth (Fig. 13).
Temperature also remains unchanged in simulation const13.The final temperatures of simulations sudden13 and lin13 are again very close, although different from 2008ss.Ice temperature changes are not restricted to the bed.In simulation sudden13, the temperature adjusts smoothly to a similar pattern as in simulation lin13; no sudden jump in basal temperature is visible in spite of the step change in basal friction parameter.Taking into account a few years for adjustment, we see good agreement for all variables after this 13year period with simulations sudden13 and lin13.This suggests that the simulations are robust towards details of how changes in the basal friction parameter occur.
We compare different temperature fields: (i) the steadystate temperature field obtained with the 1995 basal friction parameter (simulation 1995ss), (ii) the steady-state temperature field obtained with the 2008 basal friction parameter (simulation 2008ss) and (iii) the final temperature distribution of the 1995-2008 prognostic simulation with linear changing basal friction parameter (simulation lin13).Firstly, we observe that the results of lin13 and 2008ss show differences in the central, slow-flowing areas.This implies that the ice cap is not in steady state and decadal periods are not sufficient to reach a steady state.Secondly, the large-scale structure of temperatures is rather similar around Franklinbreen.During the 13-year prognostic simulation (lin13) a warming with respect to the initial conditions taken from run 1995ss of about 2-3 K at Franklinbreen's lateral margins as well as a slight warming in its fastest area (texture Fig. 11a and b) can be observed.Those 13 years of evolution are not sufficient to reproduce the cold area between the two branches visible in the steady-state simulation 2008ss (Figs. 6 and 11a).A continuation of the lin13 simulation beyond 13 years with constant friction parameter shows that a steady-state result similar to the run 2008ss is reached after about 100 years (result not shown).This is a quicker equilibration than the central area.We explain the colder areas of the steady-state 2008ss in comparison to the end of the 13-year prognostic run lin13 by the fact that the friction heating-thickness feedback is compensated by cold ice advected into the fast-flow region from the catchment area.Under the assumption that a temperate base is needed to allow for sliding, it becomes clear that another mechanism is needed to sustain the temperate base underneath the fast flow as already observed in Sect.5.3.1.
In what follows, we limit the discussion to one of the simulations 1995-2008 including firn heating (lin13firnss rather than lin13firnevol) since the transient effects are insignificant, not reaching the base within relevant timescales.The additional heat source due to firn heating does not significantly alter velocities and surface elevation or how the different variables adjust to changes in the basal friction parameter (sudden jump versus smooth adjustment).The final strain heat distribution is very similar, though friction heat is slightly greater further up in the catchment areas of the outlets (Fig. 13).The "sliding area", defined by temperatures above 271 K, i.e. temperatures close enough to pressure melting to allow for sliding (Hindmarsh and Le Meur, 2001), remains almost constant in shape and size (Fig. 11b).It is, however, challenging to interpret the changes occurring to this "sliding area", which is essential for temperature-based parameterisations of sliding laws as implemented by Seddik et al. (2012) or Dunse et al. (2011).In such laws sliding coefficients are mainly defined by the existence of a temperate base, but they are not based on physical mechanisms.They would require a precise knowledge of the temperature distribution as initial condition which is highly uncertain for ice caps like VSF.
As in the steady-state simulations including firn heating firnss, a good correlation between fast flow and temperate base is produced throughout the full 13 years.This leads us to the conclusion that firn heating is the main driver for maintaining fast flow by supplying the fast-flowing outlets with enough warm ice from the interior of the ice cap.We speculate that this convected front of increased ice temperature from firn heating had already reached the southern outlets well before 1995, explaining the observed fast velocities, but only reached Franklinbreen at around 1995.
There are several challenges to our hypothesis.First of all, we are not able to judge to what extent the initial quiescent state, given by the steady-state solution 1995ss that excludes firn heating, reflects the de facto existing englacial temperature field at this time.Secondly, a steady state including firn heating (firnss, Fig. 9) certainly overpredicts the contribution of warm ice advected to the outlets, since Franklinbreen was not yet fast flowing in 1995.
We conducted an additional simplified transient simulation starting from the initial state of 1995ss with fixed friction parameters (beta95T = those obtained for 1995), but with firn heating that on infinite timescales would lead to the same distribution as obtained with the steady-state run firnss instantly applied to the system.This simplistic simulation neglects all climatic changes leading to variations in firn heating during the last centuries as well as the heat redistribution through meltwater.The evolution of the area of warm base is depicted in Fig. 15.We observe that changes in englacial temperature caused by excess firn heating take a few centuries before a significant area of the base of Franklinbreen becomes temperate.We equally observe that different areas of the ice cap approach the steady state firnss with different timescales.This implies that it is possible that temperature at the drilling hole location could have been quite close to firnss in the 1990s, while it was not the case over the whole area or the catchment area of Franklinbreen.Typical timescales (several centuries) for this advection process to reach the lower part of the outlets can be deduced from Fig. 14.By choosing a fixed set of friction parameters, these values certainly are not exact, especially since the feedback between warming and increase of sliding velocity and hence shorter advection times is neglected.Consequently the advection times for the southern outlets are certainly underestimated since they are calculated with already high sliding velocities in these outlets.However, it provides us with a crude estimate of the advection timescales involved.These timescales (several centuries) exclude the possibility that changes in firn heating during the 20th century (as for example in firnevol) could be responsible for the trigger of fast flow at Franklinbreen.Earlier evidence of refreezing is known from the temperature profile measured at the deep borehole (steady-state firn heating run firnss).Figures 14 and 15 show that the associated heat release could have started to contribute to maintain a temperature base in the late 1990s.Consequently, we would expect that at some point the acceleration at Franklinbreen will stop and velocities level out, just as occurred decades earlier on the southern outlet glaciers of VSF.One explanation for the vast delay of Franklinbreen can be seen from Figs. 14 and 15, which document a significantly longer (in comparison to other outlet glaciers) convection timescale for ice from the interior to reach Franklinbreen.This is in line with greater distance of the firn line relative to Franklinbreen (Fig. 9).In addition, we cannot conclude about the role of a possible initial trigger for the fast flow through a thickness-friction heat feedback and resulting changes in advection times compared to the arrival of the convected front of warm ice.Furthermore, several distinguishing features of Franklinbreen have to be taken into account: Franklinbreen might have been surging in 1956 (Hagen et al., 1993), even though this was questioned by Sneed (2007) because of the lack of geomorphical evidence.It has a much bigger catchment area and probably a larger hydraulic head.It is longer and probably flatter.The flow is split into two branches and additionally complicated by a Nunatak.It terminates in a long cold fjord with more fresh water and is less prone to calving.
The question about the delayed speed-up of Franklinbreen could only be fully answered by including hydrology in the model to get a full picture of the water and heat redistribution throughout the ice cap, in particular since the englacial temperature distribution affects the englacial hydrology and vice versa.

Conclusions
We present the basal friction parameter distribution obtained on the VSF ice cap for 1995, 2008 and 2011 from the inversion of observed velocity fields using Elmer/Ice.An acceleration between 1995 and 2008 on the outlet Franklinbreen is reflected in a pronounced temporal variation of the distribution of the basal friction parameter.The drastic change of inversely determined basal friction parameters between these snapshots (Sect.5.1) renders prognostic simulations using a temporally constant parameter distribution highly inaccurate.Instead, models incorporating the full complexity of physics involved in basal sliding (hydrology, thermodynamics, tilldeformation) have to be developed for simulations of future scenarios of the ice cap.Coupling the inversion to temperature simulations affects the obtained basal friction parameter distributions only marginally because of the small influence of the viscosity-temperature dependence.
The temperature profile measured in a deep borehole in the accumulation zone has been successfully modelled and qualitatively explained by a recent change in climatic forcing at the surface influencing the firn heating (Sect.5.2).Qualitative interpretations of our firn heating experiments are limited by the uncertainty in the basal friction parameter distribution which comes into play through its influence on advection.Recent deviations from an equilibrium profile triggered by increased firn heating can explain the current shape of the upper part of the profile, whereas the much earlier existence of firn heating explains the warm temperatures observed deeper in the ice in the central thick slow-flowing areas of the ice cap.Simulated timescales for temperature changes caused by increased firn heating due to changing surface conditions to reach the ice base are upwards of a century.
A temperate base is a prerequisite to allow for sliding in the fast-flow areas.Friction heating and strain heating significantly contribute to heat production in the fast-flowing outlets, the former superseding the latter.Nevertheless, in steady-state simulations that account for only these two heat sources, Franklinbreen was the only outlet to develop a temperate base (Sect.5.3.1).Taking observed surface elevation changes into account, we conclude that a combined thickness-friction heating feedback could be the trigger for the onset of Franklinbreen's recent acceleration.In order to reach pressure melting point at the base of all outlet glaciers, we need to include the effects of firn heating in the accumulation zone and the subsequent transport of the heat into the fast flowing outlets by advection (and possibly meltwater).Thus, we identify firn heating as the heat source responsible for the persistence of the observed fast flow over all outlets (Sect.5.3.1).Prognostic simulations (Sect.5.3.2) confirm that the inflow of such warm ice is necessary to maintain a temperate base.From these simulations we also conclude that the englacial and basal temperature distribution is not in equilibrium; timescales are estimated at decades to centuries to reach steady state.However, because of the low contribution to velocity of ice deformation compared to sliding, the impact of the temperature dependence of the viscosity is small.This justifies the thermal steady-state assumption during the inversion for the sliding parameter or more generally for the initialisation of temperature fields in prognostic simulations, where the focus is not on an accurate and detailed temperature regime.Various speculations on the delayed onset of the fast flow of Franklinbreen compared to the southern outlets can be made.This question could only be fully answered when including hydrology in the model to get a full picture of the water and heat redistribution throughout the ice cap.

Figure 1 .
Figure 1.Surface velocities from remote sensing data in December/January 1995/96, December 2008 and December 2011 (original data, left side).To the right, the location of the ice cap in Svalbard (middle), surface (top) and bedrock topography (bottom) are illustrated.Sea level is indicated by a bold line.FB indicates the location of Franklinbreen, SB Sabinenbreen, RB Rijpbreen, BB Bodleybreen, AB Aldousbreen, FzB Frazerbreen, IB Idunbreen and GB Gimblebreen.
is based on the minimisation of a cost function when solving the Stokes equations iteratively with two different sets of boundary conditions.The definition of the cost function and the minimisation algorithm follow Gillet-Chaulet et al. (2012) and Jay-Allemand et al. (2011).

Figure 2 .
Figure 2. Distribution of the basal friction parameter in 1995 when coupling to temperature: (a) simulation beta95T.The other three figures show a close-up over Franklinbreen (using the depth-dependent temperature profile) in 1995 (b), 2008 (c) and 2011 (d) (simulations beta95, beta08, beta11).Note that only the southwestern corner has been patched with 1995 data in the 2008 and 2011 data sets.All figures are in logarithmic scale (units are MPa yr m −1 ).

Figure 3 .
Figure3.Evolution of the total cost function, J tot -Eq.(13), log scale, scaled by 10 9 -when iterating inverse method and temperature steady-state calculation (1995 velocities, simulation beta95T).The graph is cut at 1.6 × 10 9 for better visibility, higher values in the first iterations are thus not visible.In orange (void circles), the evolution of the cost function when not coupling with temperature is shown.

Figure 4 .
Figure 4. Evolution of maximum (and minimum) values of cumulated mesh adjustments and vertical velocities during the surface relaxation surfrelax95 to determine the end of the procedure.

Figure 5 .
Figure 5. Vertical corrections (in metres) of the topographic data (1995 basal friction parameter field) at the end of the surface relaxation (simulation surfrelax95): (a) absolute corrections, (b) relative to the initial ice thickness.

Figure 6 .
Figure 6.Basal temperature field ( • C, relative to pressure melting): (a) advection only (1995ssA), (b) adding of strain heating (1995ssSH), (c) adding of friction heating for 1995 (1995ss), (d) advection, strain and friction heating for 2008 (2008ss).The position of the ice core is indicated in red, the position of the two locations used in Fig. 7 in white.

Figure 8 .
Figure8.Temperature profile as measured (red) in a drill hole(Motoyama et al., 2008) (temperatures relative to pressure melting), the position is indicated in Fig.6(red dot).The dotted black line corresponds to the depth-dependent profile, the green line corresponds to the modelled temperature field including only strain and friction heating.The dark blue line is the result of the equilibrium firn heating simulation with the first set of parameters (percolation depth of 13.2 m, simulation firnss), which acted as initial profile for the succeeding transient runs.The orange, cyan and pink profiles illustrate the evolution of temperature in the succeeding time-evolving simulation with the second set of parameters (percolation depth of 17.6 m) evolving towards a new equilibrium (grey line).The best fit to the data among these profiles after 35 yr is highlighted in black (simulation firnevol).

Figure 9 .
Figure 9. Basal temperature field ( • C, relative to pressure melting,1995): adding of firn heating (simulation firnss).This figure remains unchanged whether the non-equilibrium part (simulation firnevol) of the firn heating is added or not.The position of the ice core is indicated in red, the position of the two locations used in Fig.7in white.The firn line is drawn in purple, the 20 m ice thickness isoline in light grey.The latter is a good indicator for underestimated ice thickness in areas where fast sliding seems to be in strong contradiction with too cold basal temperatures.

Figure 10 .
Figure 10.Change in surface elevation, S(m), during the prognostic 1995-2008 simulations for the three scenarios: (a) constant basal friction parameter, simulation const13, (b) sudden jump after 5 years, simulation sudden13 and (c) linear change in basal friction parameter, simulation lin13.

Figure 11 .
Figure 11.In a close-up of the Franklinbreen area the change in the extent of the "sliding area" (basal temperatures above 271 K) is shown.In (a) the texture is the 1995ss temperature distribution (temperatures relative to pressure melting).The light blue corresponds to the extent of the "sliding area" for 1995ss, dark blue for 2008ss and pink for the end of lin13.In (b) the temperature at the end of the 13 yr simulation when neglecting firn heating (lin13) (temperatures relative to pressure melting) is used as texture -same colour scale as (a).The light colours present the "sliding area" when neglecting firn heating (lin13), dark colours when including it (lin13firnss); blue lines represent the start of the 13 yr simulation and pink the end.

Figure 12 .
Figure 12.(a) Difference between surface and bedrock velocities, i.e. the velocity due to deformation in 1995.Vertically integrated strain heat (b) and friction heat at the bedrock (c) in the simulation 1995ss.

Figure 13 .
Figure 13.Close-up over Franklinbreen.Distribution of the mechanical heat sources at the end of the prognostic simulation.The upper line shows strain heat, the lower friction heat.Results of the simulation without firn heating (lin13) are shown in (a) and (c), the distributions are identical to those of 2008ss.Results with firn heating (lin13firnss) are shown in (b) and (d).In pink the contour of 271K and in yellow the iso-velocity line of 150 m yr −1 is drawn.

Figure 14 .
Figure14.Typical advection timescales from the region above the firn line in the centre of the ice cap to the lower parts of the outlets estimated from backward streamlines obtained with the backward Runge-Kutta method implemented in ParaView(Ahrens et al., 2005) assuming the basal friction parameter beta95T.The trajectories are coloured accordingly to the advection time to its lower end (blue colour).The firn line is depicted in orange, as texture the temperature field 1995ss is drawn.

Figure 15 .
Figure 15.Temporal evolution of the temperature field during a simplified transient simulation from 1995ss towards the steady-state firnss with the friction parameter beta95T.The firn line is shown in violet, the drill hole corresponds to the red dot.The evolution of the "sliding area" (271 K) from 1995ss to firnss is shown in rainbow colours (black area; blue, green, yellow, red contour lines) 100 years apart (dark green for example corresponds to 500 years).As texture the surface elevation is depicted.

Table 1 .
Numerical parameters used in this study.

Table 2 .
Summary of simulations.