The Cryosphere Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lov́enbreen , Svalbard

We present steady state (diagnostic) and transient (prognostic) simulations of Midtre Lov énbreen, Svalbard performed with the thermo-mechanically coupled full-Stokes code Elmer. This glacier has an extensive data set of geophysical measurements available spanning several decades, that allow for constraints on model descriptions. Consistent with this data set, we included a simple model accounting for the formation of superimposed ice. Diagnostic results indicated that a dynamic adaptation of the free surface is necessary, to prevent non-physically high velocities in a region of under determined bedrock depths. Observations from ground penetrating radar of the basal thermal state agree very well with model predictions, while the dip angles of isochrones in radar data also match reasonably well with modelled isochrones, despite the numerical deficiencies of estimating ages with a steady state model. Prognostic runs for 53 years, using a constant accumulation/ablation pattern starting from the steady state solution obtained from the configuration of the 1977 DEM show that: 1 the unrealistic velocities in the under determined parts of the DEM quickly damp out; 2 the free surface evolution matches well measured elevation changes; 3 the retreat of the glacier under this scenario continues with the glacier tongue in a projection to 2030 being situated ≈500 m behind the position in 1977. Correspondence to: T. Zwinger (thomas.zwinger@csc.fi)


Introduction
Midtre Lovénbreen (78.53 • N, 12.04 • E) is an alpine-type valley glacier in northwest Spitsbergen.Like most of Svalbard's glaciers, it is polythermal.This means that the glacier has a layer of cold ice (temperature below melting point) on top of the temperate ice layer (temperature at the melting point) (e.g., Björnsson et al., 1996;Moore et al., 1999).The glacier is representative of many small glaciers in the Arctic archipelagos, and hence could be considered to be of somewhat wider interest than may be expected considering its modest size.The glacier is also very conveniently placed about 5 km from the Ny Ålesund research station and has received a great deal of attention from glaciologists over the last 20 years.These studies provide a useful collection of results that can be used to verify, modify and validate ice flow models.
Midtre Lovénbreen has a warm-based core, cold-based ice around its snout and margins and cold ice at its surface.It has a number of basins feeding the main glacier tongue.In 2000 AD the glacier was circa 6 km long and has a maximum thickness of circa 185 m.Its elevation is 50 m above sea level at the snout and about 600 ma.s.l. at the headwall.Over 35 years the equilibrium line altitude averages circa 395 ma.s.l.Annual mass balance measurements conducted by the Norsk Polarinstitutt (NPI) since 1967 (Lefauconnier et al., 1999;Kohler et al., 2007) show consistently negative balances in all but 2 years (1987 and 1994).This is consistent with the continuous retreat of the glacier terminus by a total of about 1 km from its neoglacial maximum in the early 20th century (Hagen et al., 2003).Since 1969, mass Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Zwinger and J. C. Moore: Simulations of Midtre Lovénbreen balance measurements have been made at a series of stakes along the glacier centre line (Rippin et al., 2003).As with many or all polythermal glaciers, superimposed ice plays a major role in mass and energy balance.Wright et al. (2005Wright et al. ( , 2007) ) used a surface energy model with an explicit calculation of meltwater refreezing, together with observations of superimposed ice formation (Wadham et al., 2006) to show that superimposed ice accounts for an average 37% of the total net accumulation under present conditions.Wadham et al. (2006) used a three dimensional steadystate model to investigate the distribution of superimposed ice along the glacier centre line.The model was based on the first-order algorithm of Blatter (1995) and aimed to reproduce the contemporary three-dimensional flow regime (i.e., stress/(strain) fields) and to trace a series of synthetic annual layers introduced under a variety of mass balance regimes.The model was applied to the bed and surface geometry described by Rippin et al. (2003).As with Wadham et al. (2006), we use annual surface temperature distribution based on 30 years of measurements made by Det Norske Meteorologiske Institutt at sea level in Ny Ålesund.They give a mean annual temperature of −6.0 • C at sea level.Measurements of air temperature on the glacier from 1998-2002 (Wadham and Nuttall, 2002), indicate a mean annual temperature lapse rate of 0.004 • Cm −1 .We chose the temperature distribution in the upper snowpack such that it has the same lapse rate, but with an offset of 3 • C, such that the value of about −4.5 • C measured during the 2002 melt season at stake 7 (Wright et al., 2007) is recovered.Wadham et al. (2006) used their flow model to map the temperate region of the glacier by estimating where basal sliding was important in surface velocity measurements along the centre line.The region of inferred sliding agreed well with the region of warm basal ice inferred by Björnsson et al. (1996); Rippin (2001) on the basis of multiple-frequency radio echo sounding, and is also consistent with radar sounding done in 1998 by Moore that was used to create the bedrock DEM.Elsewhere along the centre line, the modelled internal deformation component accounted wholly for the observed surface velocity regime to within ±0.25ma −1 .
In this paper we outline a simple but useful approach to dealing with the polythermal nature of the glacier and the impact of superimposed ice on the energy balance of the glacier.We then outline the creation of the mesh boundary surfaces in the form of the free surface and bedrock DEM, the surface mass balance and geothermal heat flux estimate.The model is run in both steady-state conditions and timeevolving modes.We show how the models compare with surface evolution of the data from detailed LiDAR mapping of the 2005 and 2003 surfaces elevation and velocity field (Rees and Arnold, 2007).

The influence of superimposed ice formation
The glacier is marked by 9 stakes (numbered 2-10 for historical reasons) along the centre line (Fig. 1.).These stakes represent locations of ablation measurements and have been frequently located by geodetic GPS.The glacier is characterised by three distinct regions 1. a region below stake position 9 ( 400 ma.s.l)where the surface heat balance can be characterised by an annual mean temperature distribution for the uppermost region of the ice/snowpack of the form where z is the elevation in metres and the temperature T is given in • C.This region is essentially the ablation area of the glacier where solid ice rather than firn lies below any seasonal snow pack.Wadham et al. (2006) suggest that this region may, especially in the upper regions of the glacier, contain superimposed ice more than The Cryosphere, 2. a region where superimposed ice formation governs the surface heat balance, between stake positions 9 and 10 (∼400-450 ma.s.l).This region either gains or looses mass each year depending on the relative amount superimposed ice formed each year as a result of melting and re-freezing on the glacier surface.Wadham et al. (2006) suggest that superimposed ice in this region is recent (less than 100 years old).Later we will show how a volumetric heat source accounting for the release of latent heat by refreezing is introduced in addition to the lapse rate temperature distribution as for region 1.
3. a completely temperate region above stake position 10 ( 450ma.s.l)where the glacier is at the pressure melting point throughout its whole depth except for the seasonally cold winter snowpack.
In the computational mesh, the region boundaries were inserted manually.Consequently, their positions with respect to elevation are approximated within a range of a few metres.Those different regions and the surface temperature distribution obtained from the diagnostic run are depicted in Fig. 1.Initial model runs not accounting for the additional release of latent heat caused by refreezing processes produced unrealistic bedrock temperature distributions in region 2, where the area of temperate base -in contradiction to measurements -is significantly reduced.As the change in bulk temperatures is about one • C, velocity distributions are far less altered by neglecting the additional heat release caused by refreezing, provided its contribution to the mass balance implicitly is accounted for.Wright et al. (2007) made an intensive comparison between different approaches on how to model the impact of refreezing on Midtre Lovénbreen.Their main criterion in evaluating the approaches was comparison with measured mass balance in regions of superimposed ice formation.In our study, both diagnostic runs and prognostic models with prescribed surface accumulation were run, with the mass balance being a result of the fixed geometry or the kinematic boundary condition, respectively.Hence the main focus here is in getting a physically realistic model for the heat-transfer.As the applied integration times step sizes are longer than the size needed to resolve refreezing process time-scales, we rely on time average approaches.This led us using the Wright P max -model for formulating the impact of superimposed ice.This approach utilises the different characteristic shapes of temperature profiles as a function of depth, d, before seasonal superimposed ice formation and thereafter, Equations ( 1) and (2) contain three parameters: the mean annual air temperature, T a , the mean winter air temperature, T w and the typical penetration depth of the annual temperature cycle into the glacier ice, d ice .We chose the values of T a =−4.5C • , T w =−6.3 • C and d ice =5m as they were obtained for the region around stake position 7 (Wright et al., 2007).With the heat capacity of ice, c, taken to be constant c(T a )=2097Jkg −1 K −1 and the density of ice ρ ice =910kgm −3 , we obtain the averaged annual energy release per unit volume from superimposed ice formation as a function of depth e ice = cρ ice T summer (d) − T winter (d) . (3) The volumetric heat source (in Wm −3 ) within the layer of superimposed ice formation then is given by Eq. ( 3) divided by the amount of seconds per year 2 The numerical model Numerical simulations have been made using the open source Finite Element Method (FEM) code Elmer (http: //www.csc.fi/elmer).Elmer solves the thermo-mechanically coupled problem taking all stress components into account, i.e., solving the full-Stokes problem of the flow.As FEM utilises unstructured grids, a variable mesh size to account for the resolution can be applied (e.g., higher resolution in regions of steeper surface slope).The size of the computational mesh also necessitated use of Elmer's capability for parallel computation by using domain decomposition.

Geometry and meshing
A surface DEM was created from a digitised form of the Norwegian Polar Institute map (1979) which was constructed from aerial photographs taken in 1977.The glacier bed topography was mapped by ground penetrating radar in 1998 (Rippin et al., 2003).The glacier margin was taken from NPI maps of the Brøgger Peninsula (Norwegian Polar Institute, 1990).
A digital elevation model (DEM) of the bedrock was created from radar sounding data and measured positions of the glacier outline Kohler et al. (2002) using MATLAB.A radial bias interpolation was used to produce a DEM consisting of equally distributed points over the whole area.Due to a lack of data points, the bedrock in the upper eastern part of the glacier was partly extrapolated.As will be discussed later, this led to unphysical high velocities in those parts.
Based on the surface DEM, a footprint of the glaciated area was defined.The projection of this footprint in the horizontal plane was meshed using the commercial meshing tool Gambit.A feature of Gambit allows various mesh size functions to be utilised to produce higher horizontal mesh density www.the-cryosphere. in regions where physics or geometry demand higher resolutions, such as at steep margins.Using the DEMs of the free surface and the bedrock as inputs to a surface interpolation scheme, the three-dimensional mesh was created by extruding the two-dimensional footprint from the interpolated bedrock to the free surface, producing 10 layers of varying thickness over the whole glacier surface.The resulting mesh has a total of 120 000 nodes (135 000 brick elements).Due to the size of the problem, a domain decomposition was made to allow parallel computation on four processor cores (see Fig. 2).With such a set-up, a diagnostic run only takes 1-2 h (depending on accuracy demanded) and the 53 year transient simulation finishes within one day on a modern system of two dual-core AMD Opteron processors.
Starting from the footprint mesh, the previously discussed three distinct thermo-dynamic regions were built into the mesh.Additionally, a vertically aligned flow section of the glacier that includes all surface stake positions was introduced for post processing reasons.

Field equations and their implementation
On long time scales, the flow of ice is governed by the Stokes equation where v=(u,v,w) T is the velocity vector, p the pressure, ρ the density and g is the acceleration caused by gravity.
Density, ρ, is assumed to be constant equal to that of solid ice, ρ ice =910kgm −3 .This is a particularly reasonable assumption on Midtre Lovénbreen where the firn region is very small indeed and seasonal snow packs are of order 1 m thick.The deviatoric stress tensor, τ , is related to the deviatoric part of the strain rate tensor, D, with its second invariant, δ= 2tr(D 2 ), via Glen's flow law with an exponent n = 3) The flow rate factor, A(T ), is derived from the Arrhenius law (Paterson, 1994).Two different sets of simulation runs with different enhancement factors, E=1 and E=3, were undertaken.As the latter resulted velocities one order of magnitude too high (up to 25 ma −1 ) compared with measured values, the constant value E=1 was set.
We use a no-slip condition for the velocity at the bedrock, v bed =0.At the free surface, a vanishing deviatoric stress vector, t, is imposed.Further, the atmospheric pressure and its change over the glacier surface are taken to be zero p surf.= p atm. ≈ 0, grad p surf.≈ 0.
Here, n denotes the outward pointing surface normal.Equation ( 7) is implicitly fulfilled by imposing the natural boundary condition (defined for the Cauchy stress tensor) for the variational form of Eq. (5).The system is solved using The Cryosphere, 3, 217-229, 2009 www.the-cryosphere.net/3/217/2009/first order functions and stabilized finite elements (Franca et al. , 1992) in combination with a direct parallel solver.
For the transient (prognostic) runs, an additional kinematic boundary condition for the elevation of the free surface, s, with a given accumulation/ablation function, a ⊥ , has to be solved.The variational form of Eq. ( 8) is obtained analogously to Gagliardini et al. (2006).However, in contrast with their approach the variational inequality resulting from Eq. ( 8) under the constraint where b is the elevation of the bedrock, is not solved using a penalty method.We instead use the method for solving the contact problem as introduced by Durand et al. (2009).
In order to also impose Eq. ( 9) on the domain boundaries (the outline of the glacier), the Dirichlet boundary condition s outline =b+h min was set.Here we set h min to 5 m.As a consequence of Eq. ( 9), even deglaciated areas retained a minimum flow depth of h min , which allowed the horizontal positions of the mesh over the bedrock to be kept fixed in transient runs.Too small a choice of h min would lead to highly distorted elements and hence interfere with numerical stability.On the other side, a too high value would significantly alter the geometry of the domain.Comparison of purely mechanical test runs between the standard mesh and a mesh with h min set to 10 m revealed no significant difference in the dynamics.This and the good match between measured and computed velocity fields gives us confidence that the choice of an Eulerian grid in the horizontal mesh for the free surface evolution does not have a strong impact on the solution.
In transient simulations the evolution of the free surface implies deformation of the computational mesh.Consequently, the mesh has to be vertically rearranged after each time step.This is done by solving a virtual elasticity problem, with the mesh deformation being the nominal variable.
The temperature, T , in the heat transfer equation, is constrained by the upper limit imposed by the pressure melting point Since the glacier is everywhere less than 200 m thick, we can assume a constant 1-atmosphere pressure melting point of T pm =0 • C. The material parameters, heat capacity c and heat conductivity κ, occurring in Eq. ( 10) are functions of the temperature.This, and the variational inequality resulting from Eq. ( 11), make the solution of the heat transfer equation a non-linear problem, to which we apply the same method as explained by Zwinger et al. (2007).Throughout the glacier a volumetric heat source σ =σ strain +σ ice was calculated, where strain heating, σ strain =τ :D, and σ ice is the latent heat component Eq. ( 4) which is zero except in the superimposed region (see Fig. 1).The boundary conditions for the temperature at the free surface are implemented according to the three distinctive zones discussed earlier.At the bedrock, a standard value of 56.05 mWm −2 (Björnsson et al., 1996) was taken for the geothermal heat flux, q geo.=κ grad T •n| bed .
The solution for the age of the ice, A, is obtained by solving the purely advective equation with the boundary condition at the free surface The surface normal, n, by definition points outwards.Consequently, the condition in Eq. ( 13) guarantees that zero age is only set at inflow, whereas at outflow -by the nature of the applied discontinuous Galerkin method (Brezzi et al. , 2004) -the age is provided as a solution.As our integration times are limited, we are forced to look at steady state solutions -even the prognostic runs are essentially 50 year transient runs not fully time evolving over millennial scales required to simulate full glacier-climate interactions.Thus the first contribution on the l.h.s. of Eq. ( 12) is set to zero.

Solution strategies
The initial glacier geometry was given by the DEMs with the surface pertaining to the year 1977.The inherent problem if starting from a DEM is, that an instantaneous solution for both, the mechanical and the thermo-dynamical problem is needed as a starting point.The thermal initial conditions are critical due to the high storage capacity of the ice for heat and the low advection/diffusion rates on the glacier.In other words, it is not possible to make a correct spin-up of the system, which would demand a transient run starting from deglaciated conditions with a long enough spin-up time and -in consequence -needs realistic forcing (temperature, net accumulation).
In addition, starting from a zero velocity field at a given geometry, as is the case of a diagnostic run, is a numerical challenge, especially if dealing with a full-Stokes code in combination with a shear-thinning flow law that shows a singularity at vanishing shear-rates.In order to avoid these problems and to obtain reasonable initial fields, the following sequence of simulations is conducted: 1. System initialisation: a steady state run with a temperature and shear independent (linear) viscosity, η=10 13 kgm −1 s −1 , is set up.Using the obtained velocity solution, the heat transfer equation is run in a postprocessing step to get an initial guess of the temperature distribution.
www From the age distribution patterns it is clear that from stake 2 to the terminus the stake line (white annotated line) strongly deviates from a flow line.

Thermo-mechanically coupled diagnostic simulation:
Using the results of the previous initialisation run as initial guesses, Eqs. ( 5) and ( 10) are solved in a fully thermo-mechanically coupled manner.The shear viscosity, η(δ,T ), occurring in Eq. ( 6) now accounts for non-linearity as well as its dependence on temperature.As a consequence of the assumption of steady-state, the free surface is kept constant at the values given for the 1977 DEM, s=s 1977 .

Thermo-mechanically coupled prognostic simulation:
Using the results of the previous diagnostic run as the initial condition for all unknowns, a transient run is performed.The thermo-mechanically coupled system Eqs.( 5) and (10) as well as the time-evolution of the free surface equation ( 8) is solved on each time-level.
As will be shown in the following section, the final step is usually necessary, as errors in DEMs and also non-steady effects demand a temporal adaptation of the free surface in order to yield realistic solution, especially for the mechanical part of the problem.

Diagnostic results
Starting from the initial velocity field from step 1 listed in the previous section, obtained with a linear flow law, a thermomechanically coupled simulation was run in steady state (assuming all time derivatives to vanish).The Stokes solver and the heat equation solver were run using a fixed point iteration scheme as long as the relative residual of all solution variables fell under a given threshold.
The surface velocities found are highly unrealistic (first picture in Fig. 11), and clearly come from errors in the DEM.Investigation showed that the unrealistic velocities come mainly from interpolations in an area of the bedrock from zone 2 (see Sect. 1.1, super-imposed ice formation) that was poorly covered by the radar sounding.In this region, flow velocities locally near stake 9 exceed 15 ma −1 (4.1 cmd −1 ), which is in contradiction to measured values that are less than 5 ma −1 (1.4 cmd −1 ) (e.g., Rees and Arnold, 2007;Wadham et al., 2006).This leads to both, the age (see Figs. 3 and 4) and temperature distributions (see Figs. 5 and  6), being suspect in that particular region of the glacier, since both these parameters are highly dependent on advection.The pattern of basal temperate ice in Fig. 5 shows that the western side of the centre line has a greater region of basal melting than the easterly half of the glacier, especially towards the snout of the glacier.This pattern is well-matched by the character of the basal radar reflection (Rippin, 2001).Despite the distortions caused by the errors in the DEM, clear results are, that the basal thermal condition (Figs. 5 and 6) and to a large extent the age of the ice (see Figs. 3 and 4) depends on advection from the region of temperate ice at the uppermost parts of the glacier.This interpretation is confirmed by the computed Peclet number distribution over the glacier (see Fig. 7), The Cryosphere, 3, 217-229, 2009 www.the-cryosphere.net/3/217/2009/As a cut-off value has been introduced, the gradient around stake 5 appears strongly enhanced.As shown in Fig. 3, the stake line strongly deviates from the flow-line, which explains the decreasing age towards the tongue of the glacier.where L is a typical length scale of the glacier (in our case set to L=100 m), and α=κ/(ρc) is the local thermal diffusivity of ice.The central region of the glacier surface exhibits Pe>1, suggesting advection dominates diffusion.
The reader should bear in mind that the temperature and age distribution were obtained for a steady state run.For the latter, this naturally leads to infinitely large values of the age at the bedrock, where a no-slip condition applies.If no extensive melting and drainage of ice at the bedrock allows for a steady state solution (see e.g., Zwinger et al., 2007), in the strict sense, Eqs. ( 12) and ( 13) constitute an ill-posed problem.Consequently, the solution -especially in the vicinity of the bedrock -has to be interpreted with care.In order to avoid the breakdown of the numerical model, an upper limit for the age of 10 5 years was imposed.With areas in the lower, cold based part of the glacier showing solutions at the upper numerical limit, a cut-off value had to be introduced for a better visualisation of the results.In consequence, maximum values in Figs. 3 and 4 have to be interpreted as "equal to or larger".This also leads to visual enhancement of gradients in those figures.The strong lateral variation of the age in the lower part of the glacier can be a consequence of  15), (black) as a function of distance along the stake line from the upper glacier.
ice flow from different tributaries (and consequently different age) merging in this part.A second run with a different placement of the domain boundaries lead to the same age distribution, which excludes a numerical artifact arising from domain decomposition.The dip-angle, α, of the isochrones, is displayed in Fig. 8. Internal radar reflection horizons have been established as isochrones in ice sheets and glaciers (Eisen et al., 2003;Pälli et al., 2003).In ablation regions where the isochrones dip up and intersect the surface the internal reflections have been used to calibrate simple flowline models of Antarctic blue ice areas (Sinisalo et al., 2007), however, they have not been investigated by three dimensional flow models, nor in small polythermal glaciers in Svalbard.We expect that the dip angle of isochrones is sensitive both to general (average) flow conditions which tend to determine the thickness of the layers, but also depend on varying climate controlled factors such as trends in mass balance, and temperature induced flow variations.However, the sensitivity of radar reflections to changes in ice flow is largely untested.
We used the radar survey data obtained in 1997 to map dip angles along the central stake line.Figure 9 shows the dip angles found from the radar data.The simulation in Fig. 8 show large differences in angles spatially over the glacier surface.These large spatial gradients in angles are also found with depth along the central stake line.Comparison of the observed dip angles with the modelled surface dip angles along the stake line (Fig. 9) show reasonable agreement below 1500 m along the flow line, but poor agreement in the upper glacier where measured angles are below 5 • .
There are several plausible explanations for the discrepancy between the observed angles and the modelled ones: (i) The steady state assumptions in the diagnostic model; (ii) there are only have 10 vertical layers in the model, so gradients tend to be over pronounced.Thus a coarse mesh acts as a numerical differentiator; (iii) in the upper glacier there is a fairly uniform age field (Figs. 3 and 4).Hence the dip angle from Eq. ( 15) is essentially the inverse tangent of 2 rather small numbers.So while dip angles may be a very sensitive and useful way to target climate variability effects on glacier flow, the numerical considerations of the flow model used must be addressed.It may also be useful to consider mapping dip angles laterally across the centre line considering the rapid variations in dip angle produced in near the centre line of the glacier.

Prognostic results
Starting from the previously obtained steady state solution for the 1977 DEM, a transient run integrating over a total time span of 53 years up to the year 2030 has been undertaken.We used the averaged accumulation/ablation measurements starting from 1968 at the stake line (Kohler et al., 2002) and linearly interpolated between the stakes with respect to the vertical elevation, z, only.Consequently, changes of the surface elevation resulted in changes of accumulation.The accumulation pattern for the initial geometry for the year 1977 is depicted in Fig. 10.
The fully implicit first order time stepping scheme applied a time step size of one year.In order to test sensitivity to either time step sizes or spatial resolutions, runs on the same mesh with one quarter of the time step size as well as on a coarser mesh with the same time step size have been conducted.They converged to essentially identical solutions.
As can be seen in Figs.11 and 12, the unrealistic values of the velocity in the under-determined area of the bedrock DEM from the steady state solution (white circle in first picture "diagnostic" panel) are significantly damped within the first 8 years.Thereafter, the velocity distribution in general remains the same.
The simulation was also extended to times that allow for a prognostic study of the glacier evolution under the assumption that the present day accumulation pattern is continued.Figure 13 shows the extent as well as the velocity distribution at the simulation time level corresponding to the end of year 1997.The simulated outline of glaciated area at the tongue of the glacier is shown for the years 1985,1997,2005,2015 and 2030.The retreat is slightly stronger in the first half of the simulated time span, summing to about 500 m in total within the whole period between 1977 and 2030.Rees and Arnold (2007) suggest that the snout retreated about 7 m between 2003and 2005. Rippin et al. (2003) ) suggest that the snout retreated at an average rate of 8.3 ma −1 for the 18 years between 1977 and 1995.The rates between 1985 and 2005 shown in Fig. 13 average to about 7.5 ma −1 , which is consistent with the observations.
Figure 13 shows the surface velocity distribution (in cmd −1 ) and measured velocities (Rees and Arnold, 2007) from ground survey data and semi and fully automated motion sensing from DEMs made from airborne LiDAR done in 2003 and 2005 (not shown in 13).Features for reliable velocity tracking were only found in the lower (less snow covered) part of the glacier.The general pattern found by Rees and Arnold (2007) agrees with direct measurements of velocity on a stake network on the lower glacier (Rippin, 2001).The results shown in Fig. 13 show a similar pattern of increasing velocities from the terminus towards the middle of the glacier, with values very similar to those of Rees and Arnold (2007) and Rippin (2001), that is up to about 1.5 cmd −1 in the region of the glacier shown in Fig. 13.
Vertical changes in the glacier were investigated by Rees and Arnold (2007) who mapped height changes at high resolution (Fig. 14).The glacier lowered up to 5 m near the terminus, with much of the surface lowering by 1 m.More recently Barrand et al. (2009) considered errors in estimating surface height changes using the same LiDAR datasets as a function of ground control points.Their detailed analysis produced slightly less negative mass balance results than Rees and Arnold (2007), with most differences being at the high elevation parts of the glacier.The modelled thinning rates are in excellent agreement with these observations with only the very upper reaches of the accumulation basins showing rising elevations.Rees and Arnold (2007) note that there is a discrepancy between the glacier mass balance they obtain from their LiDAR surface elevation results and the glaciological mass balance found from traditional stake measurements, with the 1977-1995 stake mass balance (around 0.35 myr −1 ) being only about 70 % of the DEM mass balance.However, here we use the glaciological mass balance as input to the model and produce a mean surface lowering between 2030-1977 as 0.175 myr −1 .So although the pat-tern of lowering mirrors that Rees and Arnold (2007) and Barrand et al. (2009) observe, with maximum values being very similar, the mean over the whole glacier is considerably less negative than the DEM (0.51 myr −1 ), but closer to the long term glaciological mass balance from 1936-1962 (0.15 myr −1 , Kohler et al., 2007).This suggests that the mass loss at the snout of the glacier is not a simple long term adjustment in glacier geometry but rather to a recent acceleration of surface melting due to warmer temperatures in recent decades.We adapt a simple time-averaged (P max ) method of accounting for the impact of superimposed ice on the thermal properties of the Midtre Lovénbreen glacier.The simple superimposed ice modelling approach to polythermal Svalbard glaciers requires only temperature profiles in winter and summer to be measured rather than simulating individual refreezing events.Artifacts in input parameters to finite element models produce potentially highly degraded output velocity fields.However, relaxation of the free surface can quickly accommodate errors in input DEMs.In our case a factor of three discrepancy between steady-state modelled and observed local velocity was reduced to negligible within only 8 time steps of a prognostic simulation.Diagnostic and especially short period prognostic runs produce excellent agreement with a wide variety of observational data sets on Midtre Lovénbreen.These include surface velocity distribution, cold-temperate internal thermal boundary, glacier retreat rate and surface thinning rate.The only input data required is the initial surface and bedrock DEMs, the geothermal heat flux (taken as a standard value), and a surface energy balance from the air temperatures and superimposed ice parametrisation.However, field knowledge of the glacier is essential to delineate the three different zones of the glacier (firn region, superimposed ice and ablation areas).The Peclet number map shows that advection dominates diffusion over the majority of the glacier surface, implying that accurate mass balance must take into account horizontal motion as well as simple vertical glacier movement.
The age distribution in the glacier is an interesting parameter as old ice is certainly available over most of the glacier surface, which may be useful for certain paleoclimate purposes.The surface age distribution appears to be rather sensitive laterally across the centre line, with the oldest ice being confined to a narrow band along the centre line.Preliminary comparisons of radar internal reflection horizons and modelled isochrones dip angles indicate encouraging agreement over much of the glacier.However, there are important numerical considerations that may produce artifacts in dip angles, particularly relatively coarse vertical resolution.However, radar dip angles are probably a very sensitive tool to examine climate-glacier dynamical links, on timescales of millennia in the case of Midtre Lovénbreen.Further development and investigation requires a broader set of radar observations, and more critically, significantly (100 to 1000 times) more computing power to run higher resolution and longer period prognostic simulations.

Fig. 1 .
Fig. 1.The three, with respect to the surface energy balance, distinct regions of the Midtre Lovénbreen.The colour texture shows the surface temperature (obtained from the diagnostic run).The contours represent the surface elevation of the 1977 DEM.The stake positions as well as elevation iso-lines are annotated.The three different regions and their boundaries (dashed lines) are indicated.

Fig. 2 .
Fig. 2. The mesh (left) and the bedrock geometry (right) obtained for Midtre Lovénbreen.The stake positions are annotated.The colour texture of the left-hand picture indicates the domain decomposition into four partitions introduced for the parallel computations.The contours in the right-hand picture indicate the bedrock elevation above sea-level.

Fig. 3 .
Fig.3.Age distribution (with respect to 1977) of the diagnostic run 2 grid layers above the bedrock (left) and the free surface (right).Note that the the graphs have different ranges of the variable.From the age distribution patterns it is clear that from stake 2 to the terminus the stake line (white annotated line) strongly deviates from a flow line.

Fig. 4 .
Fig. 4. Age distribution (with respect to 1977) of the diagnostic run along the stake line of the glacier (colour texture, cut-off value 3000 years).The contours show isochrones in steps of 100 years for 1000 years before 1977.Vertical scales are three-fold exaggerated.As a cut-off value has been introduced, the gradient around stake 5 appears strongly enhanced.As shown in Fig.3, the stake line strongly deviates from the flow-line, which explains the decreasing age towards the tongue of the glacier.

Fig. 5 .
Fig. 5. Basal temperature distribution of the diagnostic run.The cold-temperate ice transition line is indicated by a white line.

Fig. 6 .Fig. 7 .Fig. 8 .
Fig. 6.Temperature distribution of the diagnostic run along the stake line.The cold-temperate ice transition line is indicated by a white line.Vertical scales are three-fold exaggerated.

Fig. 9 .
Fig. 9. Comparison of observed radar internal layer dip angles (blue) and modelled isochrones dip angles Eq. (15), (black) as a function of distance along the stake line from the upper glacier.

Fig. 10 .
Fig. 10.Accumulation pattern computed as a function of elevation for the year 1977.Accumulation is given in ice-equivalent with an ice/water density ratio of ρ ice /ρ H 2 O =910/1000 being applied.

Fig. 11 .
Fig. 11.Evolution of the velocity (colour texture) and the ice thickness (iso lines in 25 m intervals) of a transient run starting from the steady state solution obtained for the year 1977.The position of the glacier's terminus is marked with a red dashed line.The area where drastic changes in the free surface evolution take place (errors in DEM), is marked with a white circle in the first picture of the series.

Fig. 12 .Fig. 13 .
Fig. 12. Evolution of the velocity (colour texture) and the cold-temperate ice transition line of a transient run starting from the steady state solution obtained for the year 1977.The outline of the original surface of the 1977 DEM is indicated with a red line.The area where drastic changes in the free surface evolution take place (errors in DEM), is marked with a black circle in the first picture of the series.Vertical scales are three-fold exaggerated.

Fig. 14 .
Fig. 14. Results of the simulated change in flow depth from the prognostic simulation between the years 2003 and 2005 in comparison with the measured results from Rees and Arnold (2007) (insert on upper left corner).The black area at the glacier's tongue indicates the area deglaciated between 2003 and 2005.