Simulating melt, runoff and refreezing on Nordenskïoldbreen, Svalbard, using a coupled snow and energy balance model

A distributed energy balance model is coupled to a multi-layer snow model in order to study the mass balance evolution and the impact of refreezing on the mass budget of Nordenskï oldbreen, Svalbard. The model is forced with output from the regional climate model RACMO and meteorological data from Svalbard Airport. Extensive calibration and initialisation are performed to increase the model accuracy. For the period 1989–2010, we find a mean net mass balance of−0.39 m w.e. a−1. Refreezing contributes on average 0.27 m w.e. a−1 to the mass budget and is most pronounced in the accumulation zone. The simulated mass balance, radiative fluxes and subsurface profiles are validated against observations and are generally in good agreement. Climate sensitivity experiments reveal a non-linear, seasonally dependent response of the mass balance, refreezing and runoff to changes in temperature and precipitation. It is shown that including seasonality in climate change, with less pronounced summer warming, reduces the sensitivity of the mass balance and equilibrium line altitude (ELA) estimates in a future climate. The amount of refreezing is shown to be rather insensitive to changes in climate.


Introduction
The mass balance has been acknowledged as the critical link between glaciers and climate (Meier, 1965). Numerical models capable of simulating the surface mass balance have proved useful in analysing the temporal evolution and spatial distribution of the mass budget of ice masses in more detail than provided by observations only (Greuell, 1992;Hock, 1999). Additionally, mass balance models have been used as predictive tools to study the sensitivity of glaciers to climate change Braithwaite and Zhang, 1999;De Woul and Hock, 2005;Raper and Braithwaite, 2006;Hock et al., 2007). Distributed mass balance models that solve the energy balance to compute melt have been shown to be capable of simulating spatial melt patterns and their variability (Arnold et al., 1996;Klok and Oerlemans, 2002;Hock and Holmgren, 2005).
Refreezing of percolating and stored water in snow and firn contributes significantly to the mass balance of glaciers (Schytt, 1949;Koerner, 1970) and has a pronounced impact on the thermal structure (Greuell and Oerlemans, 1989;Jania et al., 1996). Refreezing of percolating water is most pronounced in spring when low subsurface temperatures and the presence of a snow pack increase the potential for refreezing. During the melt season superimposed ice may form when percolating water accumulates a slush layer on top of the cold impermeable ice layer and refreezes (Wadham and Nuttall, 2002;Obleitner and Lehning, 2004;Wright et al., 2005). After the melt season, water trapped in pore spaces is subject to refreezing when a cold wave penetrates into the snow/firn pack (Pfeffer et al., 1991). The latent heat release after refreezing raises subsurface temperatures and hence affects the heat flux from the surface into the ice. Refreezing below the previous year's summer surface in the accumulation zone, referred to as internal accumulation, has received considerable attention, since this term is disregarded by traditional mass balance observations (Trabant and Mayo, 1985;Schneider and Jansson, 2004;Reijmer and Hock, 2008). Other studies show the significance of refreezing in the timing and rate of englacial water transport (Pfeffer et al., 1991;Fountain, 1996;Jansson et al., 2003), which has substantial implications for basal dynamics (Zwally et al., 2002;Van de Wal et al., 2008;Schoof, 2010). Coupling of a distributed energy balance model to a snow model, which simulates the subsurface temperature, density and water evolution, is required to accurately simulate ice melt, refreezing and runoff, and to study the relative impact of refreezing on the mass budget of a glacier.
In this study, a distributed energy balance model, developed along the lines presented by Klok and Oerlemans (2002), is coupled to a multi-layer snow model based on the SOMARS model (Simulation Of glacier surface Mass balance And Related Subsurface processes) described by Greuell and Konzelmann (1994). This is then used to simulate the spatial distribution and temporal evolution of the mass balance on Nordenskiöldbreen, Svalbard. Gridded meteorological input fields are constructed from output of the Regional Atmospheric Climate Model (RACMO) (Van Meijgaard et al., 2008) and weather station data from Svalbard Airport. Extensive calibration is performed to constrain values of model parameters, whereas comprehensive initialisation is done to attain subsurface profiles at the start of the simulation. The model is run over the period 1989-2010. We present the evolution of the mass balance and energy balance and discuss this in connection with evolving subsurface properties. Sensitivity experiments are performed to investigate the sensitivity of the model output to the parameter setup, initial subsurface conditions and climate perturbations.
Mean surface temperatures at Svalbard Airport have risen by 0.22 • C per decade since 1912, and regional climate model projections for the 21st century predict a warming from  to 2071-2100 ranging from 3 • C in the southwest of Svalbard to 8 • C in the northeast . Temperature increase is expected to be less pronounced in summer than during the winter season. Regional differences and seasonal variations in climate change both have a significant impact on the sensitivity of the mass balance and should be accounted for in future mass balance projections. Sensitivity tests with seasonally dependent climate variations, based on downscaled future general circulation model (GCM) estimates, are done to investigate the relevance of accounting for seasonal variability of climate change in mass balance projections.

Nordenskiöldbreen
Nordenskiöldbreen is a large outlet glacier, situated in central Spitsbergen and connected to a large ice plateau, Lomonosovfonna (Fig. 1). Ice flows around two rock formations, Terrierfjellet and Ferrierfjellet, towards the Adolfbukta fjord. Presently, the glacier front can be divided into an actively calving part and a part which has retreated on land. The glacier front was known to be calving along its full width during at least a major part of the 20th century (Plassen et al., 2004). Calving of the retreating front may become more sig- nificant again in the next few decades, as recent Ground Penetrating Radar (GPR) observations performed during fieldwork in spring 2010 indicate bedrock heights below sea-level upstream of the current glacier front.
A Digital Elevation Model (DEM) of the glacier and its surroundings with a 40 m spatial resolution was derived from stereoscopic optical images gathered in 2007 and provided by the SPIRIT project: SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies (Korona et al., 2009). Some processing of the DEM was done to remove erroneous spikes in the surface profile. The resulting gaps were filled using interpolation techniques. A contour plot of the surface topography is shown in Fig. 1.
The area of Nordenskiöldbreen chosen, bounded by the black line in Fig. 1, extends up to the estimated position of the ice divide, where horizontal ice flow is likely to be small. It can therefore be assumed that the total mass budget of the selected grid is merely the sum of the surface mass balance and a negative mass flux by calving. The selected grid covers a total area of 193 km 2 and spans an altitudinal range of 0 to 1195 m a.s.l. The highest point of the Lomonosovfonna ice cap is located at 1237 m a.s.l. The highest ice velocities are found along the main flow line between Terrierfjellet and De Geerfjellet, where mean surface velocities are typically in the range of 30-60 m a −1 (Den Ouden et al., 2010). GPR observations indicate an ice thickness of more than 600 m at some locations along the flow line. The mean surface temperature and amount of precipitation measured at Svalbard Airport (27 m a.s.l.), situated ∼55 km west of the glacier snout, over the period 1980-2010 are equal to −6.7 • C and 190 mm per year, respectively. Owing to the high latitude, seasonality in the temperature cycle is strong, with monthly mean extrema at Svalbard Airport ranging from −16.7 • C in February to +5.9 • C in July.

Data
In this study, data are used to: (1) construct meteorological input to force the surface energy balance model (Sect. 3.1), (2) calibrate poorly-constrained model parameters (Sect. 4.4), and (3) validate model results (Sect. 5.3). Continuous measurements with an Automatic Weather Station (AWS) at an altitude of ∼524 m a.s.l. on the glacier (Fig. 1), operated by the Institute for Marine and Atmospheric research Utrecht (IMAU) since March 2009, are employed for both calibration and validation purposes.

Meteorological input
Two sources of meteorological data are used to construct gridded spatial patterns of air temperature, humidity, precipitation, cloud cover and air pressure: (1) output of the Regional Atmospheric Climate Model (RACMO) and (2) meteorological data from Svalbard Airport (SA).
Cloud cover and precipitation estimates are constructed from meteorological time-series at Svalbard Airport. Cloud cover observations with a 6-hourly resolution are downscaled to the 3-h model resolution by interpolation, whereas 3hourly values of the precipitation rate are constructed by homogeneous distribution of observed 12-hourly precipitation totals in time. On the grid cloud cover is assumed to be spatially invariant, whereas precipitation increases linearly with height at a calibrated rate of 370 mm per km (see Sect. 4.4). The precipitation rate at 27 m a.s.l. is set equal to the precipitation rate at Svalbard Airport, and above an altitude of 971 m a.s.l. the precipitation rate is assumed to be constant. The altitude of 971 m a.s.l. is chosen such that the parameterised mean maximum precipitation rate is equal to the observed mean maximum precipitation of 540 mm per year found by Pälli et al. (2002) on Nordenskiöldbreen for the period 1963-1999. Gridded 3-hourly air temperature, pressure and specific humidity input is constructed from output at a 11 km resolution of the Regional Atmospheric Climate Model (RACMO), as presented by Ettema et al. (2010). RACMO is forced at the boundaries with ERA-Interim reanalysis data, provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) for the period 1989-2010. We use data from two RACMO grid points located within 2 km from the grid at altitudes of 461 and 957 m a.s.l. Every model time-step, altitudinal gradients of air temperature, specific humidity and potential temperature between the two RACMO grid points are computed and used to linearly inter-and extrapolate these variables onto the computational grid (40 m resolution). Note that air pressure and temperature at the two RACMO grid points are used to compute the potential temperature at these points. Air pressure on the grid is then computed using gridded fields of air temperature and potential temperature. The sensible and latent heat flux formulations in the surface energy budget specifically require temperature and specific humidity input that is unaffected by the glacier's microclimate (Oerlemans and Grisogono, 2002). Therefore, RACMO humidity and temperature estimates at 32-m altitude above the glacier surface are used in the calculation of turbulent fluxes. Additionally, the 32-m temperatures are adopted in the computation of the incoming longwave radiation component. A comparison of 3-hourly gridded air temperatures at two altitudes and observed temperatures at ∼4 m above the surface at the location of the AWS is shown in Fig. 2. The observations since 2009 unfortunately do not cover the melt season due to data logger problems. High correlations are found between the observed values and RACMO 32-m (R = 0.96) and 2-m (R = 0.94) temperatures, which demonstrates that temperature variations are well replicated in the RACMO dataset. The high correlations indicate the potential to use regional climate model data to construct input fields in regions where observations are scarce. Since we lack temperature observations unaffected by the glacier surface (typically at heights >10 m above the surface), no bias correction to the employed 32-m temperatures is applied. In the RACMO domain, Svalbard is located relatively close to the boundary, which leads to significant artifacts in cloud cover and precipitation estimates in this region. We therefore decided not to use RACMO data for those variables for further analysis and use the aforementioned Svalbard Airport data instead.

Stake measurements and snow profiles
Since 2006, stake measurements have been made at multiple sites on the glacier (S1-11 in Fig. 1). Generally, stake readings are done once a year in early spring and provide estimates of surface height variations. Furthermore, snow depth measurements at the stake locations are done to estimate the relative contribution of ice melt to the change in surface height. Snow pits were dug at several sites in 2008 and 2009 in order to measure vertical profiles of density and temperature. The mean observed snow density of 372 kg m −3 is used in combination with the change in snow depth to convert the observed surface heights into mass balance estimates. Continuous surface height measurements with a Sonic Ranger in 2007 are used for model calibration, as described in Sect. 4.4. In May 1997, a 121-m long ice core was drilled near the summit of the Lomonosovfonna (DS in Fig. 1). A 36-m deep vertical density profile was constructed and presented by Pohjola et al. (2002). Borehole temperatures were measured down to the bottom of the ice core, and have been analysed by Van de Wal et al. (2002). The role of these vertical profiles in the derivation of initial subsurface conditions at the start of the simulation is discussed in Sect. 4.5.

Model setup
In this study, a distributed energy balance model is coupled to a snow model and applied to Nordenskiöldbreen. The two models are coupled in the sense that melt water production at the surface serves as input for the snow model, which simulates storage and refreezing of percolating water. Furthermore, the models are coupled through the subsurface heat flux, which affects the surface energy budget and depends on vertical profiles of temperature and density.
The specific mass balance of a certain grid cell on the glacier is defined as the accumulated exchange of mass per unit area over a period of time. It is the sum of accumulation by precipitation and riming, and ablation by runoff and sublimation. Ice melt only influences the mass budget if the produced melt water runs off and does not refreeze in the underlying snow or firn pack. Refreezing of rain water below the surface provides an additional contribution to the glacier's mass budget. Accurate computation of the mass balance of a glacier therefore involves treatment of both surface and subsurface mass fluxes.

Surface energy balance model
Forced by meteorological input data, the energy balance model calculates all energy fluxes that contribute to the surface energy budget. The sum of all fluxes is equal to the energy available for melting (Q melt ), as described by the energy balance equation: where SW net is the net shortwave radiation, LW net is the net longwave radiation, Q sens and Q lat are the turbulent sensible and latent heat flux, Q rain is the heat transfer by rainfall, and Q sub is the heat flux into the ice. Fluxes towards the surface are defined as positive. In the model, energy fluxes are formulated such that the surface temperature is the only unknown, which is found by iteratively solving Eq. (1) with the left-hand-side set to zero. In case the computed surface temperature is above melting point, the surface temperature is set to melting point and energy fluxes are recomputed. In that case the sum of fluxes is positive and melting will occur. Next, a brief description of fluxes in the energy balance model will be given. Since the model has been developed along the lines presented by Klok and Oerlemans (2002), the reader is referred to that study for further details. The amount of solar radiation impinging on the surface depends on multiple factors: the top-of-atmosphere radiation, the transmissivity of the atmosphere, shading by the surrounding terrain, the slope aspect and gradient of the grid cell, the reflectivity of the surrounding terrain and the diffusivity of the sky. Input of cloud cover, air pressure and specific humidity is required to explicitly compute the atmospheric transmissivity due to Rayleigh scattering and gaseous absorption (Kondratyev, 1969), water vapor absorption (Mc-Donald, 1960) and attenuation by aerosols (Houghton, 1954) and clouds. The model distinguishes between direct, diffuse and reflected solar radiation coming from the surrounding terrain, and determines every half-hourly time-step whether a grid cell is shaded by the surrounding topography. A halfhourly time-step is used to improve the accuracy in the computation of direct solar radiation, and values are averaged to obtain 3-hourly estimates. The computed amount of outgoing solar radiation is controlled by the formulation of the surface albedo. A parameterisation developed by  has been adopted, in which the albedo is expressed as a function of the time since the last snowfall event and the snow depth. The calculated amount of incoming longwave radiation from the sky is expressed as a function of cloud cover, air temperature and specific humidity . The Stefan-Boltzmann law, describing thermal emittance of a blackbody, is used to compute the The Cryosphere, 6, 641-659, 2012 www.the-cryosphere.net/6/641/2012/ outgoing longwave radiation. The turbulent sensible and latent heat flux formulations depend on the large-scale temperature and humidity, following the expressions given by Oerlemans and Grisogono (2002). In contrast to the Monin-Obukhov similarity theory, in this approach no knowledge is required of near-surface temperature, humidity and wind characteristics; quantities which are strongly influenced by the ice surface. Despite its small impact on the energy budget, the heat supplied by rain at the ice surface is also considered in the surface heat budget. Finally, the glacier heat flux, which depends on the conductivity of the medium and the vertical temperature gradient, is computed by extrapolation of subsurface heat transport in the two uppermost firn layers to the surface.

Subsurface model
The evolution of vertical profiles of temperature, density and water content is simulated with a subsurface model, which is based on the SOMARS model developed by Greuell and Konzelmann (1994). SOMARS has been coupled to a distributed energy balance model before by Bougamont et al. (2005) and Reijmer and Hock (2008). The subsurface temperature evolution is described by the thermodynamic equation: where ρ is the layer density, T is the layer temperature, c p (T ) is the heat capacity of snow/ice, z is the vertical coordinate, κ(ρ) is the effective conductivity, F is the refreezing rate and L M is the latent heat of melting (3.34 × 10 5 J kg −1 ). z is the layer thickness. Expressions for κ(ρ) and c p (T ) are taken from Sturm et al. (1997) and Yen (1981), respectively: The first term on the right hand side of Eq.
(2) represents heating by a vertical diffusive heat flux, whereas the second term describes heat production by refreezing of water within a layer. Since penetration of shortwave radiation in the upper snow layers is neglected, no subsurface melting can occur. This implies that all water in the snow pack originates from percolation of surface melt water or rain water. The densification equation describes the time-evolution of the vertical density profile: Here, K g (ρ, T ) represents gravitational densification, computed using a formulation developed by Arthern et al. (2010) based on in situ measurements of Antarctic snow compaction, and recently modified by Ligtenberg et al. (2011): where b is the accumulation rate (in mm a −1 ), g is the gravitational acceleration (9.81 m s −2 ), ρ ice is the density of ice (917 kg m −3 ), R is the universal gas constant (8.314 J mol −1 K −1 ), and E c (=60 kJ mol −1 ) and E g (=42.4 kJ mol −1 ) are the activation energies associated with creep by lattice diffusion and grain growth, respectively. T avg denotes the temporal mean subsurface temperature and is computed every model time-step by taking the mean subsurface temperature of the preceding year. Ligtenberg et al. (2011) introduced the dependence of C on the accumulation rate b, which yields the following formulations of C(b): In case a snow/firn pack is present, the available amount of water at the surface, originating from ice melt and rainfall, will percolate downwards. Refreezing of the percolating water will raise subsurface temperatures and densities. The density cannot exceed the density of ice, whereas the temperature cannot be raised above melting point, which exposes two constraints on the maximum amount of refreezing. In case not all the available percolating water refreezes in a layer, a small amount of water, called irreducible water, will be held by capillary forces while the remaining water percolates into the next layer. In accordance with Schneider and Jansson (2004), an empirical relation has been used to compute the maximum irreducible water content of a layer θ mi , i.e. the ratio of the mass of irreducible water to the total mass of the layer. θ mi is expressed as a function of the porosity n, i.e. the ratio of pore space to the total volume of the snow layer, as follows: Downward percolation of water continues until an ice layer is reached. On top of the impermeable ice, water may be stored in the remaining pore space to form a slush layer. In contrast to the irreducible water content, the slush water runs off gradually, resulting in an exponential decay of the slush water content S: where t is the model time-step. The time-scale t * controls the efficiency of the runoff and is expressed as a function of the surface slope β as follows (Zuo and Oerlemans, 1996): where C 1 , C 2 and C 3 are constants for which values of respectively 0.5, 200 and 133 are similar to those found by Reijmer and Hock (2008). Consequently, the runoff timescale www.the-cryosphere.net/6/641/2012/ The Cryosphere, 6, 641-659, 2012 takes values of 200, 20 and 0.5 days for surface slopes of 0, 1 and 5 degrees, respectively. Surface runoff occurs when either bare ice is exposed at the surface or when the slush water level in the firn pack equals the surface height. Under both circumstances, excess water is assumed to run off instantly. Note that a snow pack effectively delays the runoff of melt water produced at the surface, and hence the vertical transport of water to the base of the ice. Horizontal transport from one grid cell to another is not considered in the model, since this requires explicit treatment of surface flows and vertical drainage through moulins, which is beyond the scope of this study.

Numerical setup
Numerical experiments are performed on a grid containing 825 × 646 grid points with a horizontal grid spacing of 40 m. Of these points only 23 %, corresponding to 120 671 grid cells, are assigned to the glacier ( Fig. 1). Prior to the start of the time-loop, in which the surface energy balance and subsurface profiles are updated with a 3-h time-step, topographical parameters are determined following a procedure described by Dozier and Frew (1990). These topographical parameters include the terrain view factor, i.e. the fraction of the overlying hemisphere covered by terrain, and the ice fraction, i.e. the fraction of the surrounding visible terrain covered by ice. Both parameters are required in the calculation of incoming solar radiation coming from the surrounding terrain and depend strongly on the orientation and location of the grid cell. Topographic shading is time-dependent and evaluated with a half-hour time-step. Note that the glacier geometry is assumed to be invariant in time. Incorporating an evolving geometry would require coupling of the mass balance model to an ice dynamical model. By not considering the geometric evolution, we, among other effects, disregard the impact of the mass balance-height feedback, which enhances surface height fluctuations and may become significant over longer periods of time.
The snow model contains a total of 17 vertical layers extending 47 m below the surface, with layer depths increasing exponentially from 0.10 m just below the surface to about 10 m for the lowermost layer. After a change in the surface height due to snowfall, gravitational densification, ice melt or sublimation/riming, vertical profiles will shift and layer properties are updated accordingly. Diffusion of heat in the snow pack is computed by applying an explicit method, in which forward differencing for the time derivative and second-order central differencing for the space derivative is used. A sufficiently small time-step of one hour is used to assure stability and convergence of the heat diffusion method in the thin uppermost layers. At the lower boundary of the snow model, no heat transport can occur and the vertical density gradient is set to zero. Fresh snow is added at the top of the snow model with a density ρ fs of 300 kg per m 3 .

Calibration
The set of model formulations contains several poorlyconstrained parameters. Observations on the glacier provide valuable information that can help to reduce the uncertainty in these parameters. In this study, radiative flux observations at the AWS site, Sonic Ranger measurements, snow depth at two stake sites and precipitation data from Svalbard Airport are used to calibrate the model.
In order to calibrate the modelled incoming solar radiation (SW in ), AWS measurements in 2009 and 2010 are employed. Observed half-hourly SW in is corrected for the tilt angle of the sensor and averaged to obtain 3-hourly estimates. Two steps are taken to tune the modelled SW in : (1) calibration of the aerosol transmissivity under clear-sky conditions, and (2) finding an expression for the cloud transmissivity as a function of the cloud cover. The aerosol transmissivity, τ a , can be expressed as a function of the optical air mass m (pressure corrected), following Houghton (1954): where k is a constant, for which values are typically in the range 0.87-1.00, depending on the geographic position (Davies and McKay, 1989). We found a value for k of 0.97 by matching the observed SW in under clear-sky conditions to simulated values over the calibration period. The cloud fraction at the AWS site is estimated using a method described by Van den Broeke et al. (2004 and Kuipers Munneke et al. (2008). This method assumes that at a certain air temperature, LW in is linearly related to the cloud cover. Following Klok and Oerlemans (2002), the cloud transmissivity τ cl is computed by taking the ratio of the observed SW in and the modelled clear-sky SW in . The resulting values of τ cl are assumed to be a function of the cloud fraction, and similar to the procedure described by Greuell et al. (1997), a fit has been made to find the following expression for the cloud transmissivity as a function of the cloud fraction n: Comparing Eq. (12) to the expression found by Greuell et al. (1997) for an Alpine glacier shows that τ cl in Eq. (12) depends less strongly on the cloud fraction, which might be related to geographical differences in the prevailing cloud type and associated optical depth. The incoming longwave radiation flux depends on the emissivity of the sky , which is formulated as a function of the cloud fraction after Konzelmann et al. (1994): with cl a constant and cs a function of the vapour pressure e a and air temperature T a : The Cryosphere, 6, 641-659, 2012 www.the-cryosphere.net/6/641/2012/ where b is a constant. In this study, LW in is calibrated by finding values for the constants b and cl , for which discrepancies in the modelled LW in are minimized in respect to observed values. Due to data logger problems during the melt season in 2009 and 2010, only longwave radiation measurements outside the melt season are available for comparison. Note that cl and b can be determined independently by considering cloud-free and fully cloudy conditions. A value for cl was found by comparing modelled and observed overcast (cloud fraction >95 %) LW in and resulted in a value for cl of 0.990. A comparison of simulated and observed LW in under clear-sky conditions (cloud fraction <5 %) led to a value for b of 0.447. These values agree well with estimated values for cl and b of 0.984 and 0.433, respectively, by Klok and Oerlemans (2002) for the Morteratsch glacier in Switzerland.
The observed albedo at the AWS site, derived by taking the ratio of reflected and incoming shortwave radiation, is employed to constrain values for the ice albedo α ice (0.32), fresh snow albedo α fs (0.87) and to estimate a threshold snowfall rate P fs (0.04 mm w.e. per hour) above which the albedo is set to the albedo of fresh snow.
The relative amount of liquid and solid precipitation in the model is determined by the local air temperature. A threshold air temperature T s/r , for which the precipitation is assumed to be half rain and half snow, is derived from time-series of precipitation type and the corresponding air temperature at Svalbard Airport since the year 2000. The resulting value for T s/r is 274.6 K. Around this threshold value the relative fraction of snowfall with respect to the total precipitation is assumed to decrease linearly from 100 % at T s/r − 1 K to 0 % at T s/r + 1 K. Snow depth measurements at the sonic ranger site and the AWS site in March 2007 and 2008 are used to tune the altitudinal precipitation gradient, γ p . This resulted in a yearly mean value for γ p of 370 mm per km. With this gradient, the average maximum precipitation of 540 mm w.e. per year (Pälli et al., 2002) is reached at 971 m a.s.l.
The turbulent flux formulations, adopted from Oerlemans and Grisogono (2002), include a background exchange coefficient C b , associated with turbulence generated by the largescale wind. Lacking direct observations of turbulent fluxes, a value for C b of 0.0028 has been estimated by matching simulated snow and ice melt in 2007 to observations at the Sonic Ranger site. Note that calibration of C b and γ p should be done after calibration of the other tuning parameters since ice melt and snow depth do not vary independently of the setup of the other tuning parameters. We chose to tune γ p before calibrating C b , since the influence of C b on snow depth is much less pronounced than the impact of γ p on ice melt.

Initialisation
Subsurface profiles of density, temperature and water depend strongly on the history of the subsurface processes and the surface forcing. Due to the significance of vertical advec-tion and diffusion in mass and heat transport, the response time of subsurface variables increases rapidly with distance from the source of variability (mainly at the surface). At a depth equivalent to the maximum depth in the snow model (47 m), the response time is typically on the order of a few decades. It is therefore important to have properly initialised subsurface profiles at the start of the model run in 1989. To do so, observed snapshots of the vertical borehole temperatures at the drill site in May 1997 are used as a benchmark against which the initial profiles are calibrated. The strategy of the initialisation procedure is to loop the model over the period 1989-1997 as often as is required to resemble the observed borehole temperature profile. Note that the drill site is not part of the selected grid of Nordenskiöldbreen, where the surface energy budget and vertical profiles are computed every model time-step (Fig. 1). Nevertheless, for initialisation purposes the surface energy balance and evolution of the snow pack has been simulated for the grid point nearest to the drill site as well. At the start of the first loop of initialisation, linear vertical density profiles are prescribed, whereas vertical temperatures are set equal to an arbitrary value of −8 • C. The initial snow depth is set to increase linearly with altitude from 0 m at 700 m a.s.l. to 25 m at the altitude of the drill site (1237 m a.s.l.). Figure 3 shows subsurface densities and temperatures at the location of the drill site after every 8 yr loop of the initialisation procedure. Observed subsurface temperatures  and densities  at the drill site in May 1997 are also shown. Best agreement between modelled and observed subsurface temperatures is found after 32 yr of initialisation. The initial vertical profiles at the start of the fourth initialisation loop are used as reference profiles for the entire grid in May 1989. Discrepancies in the simulated temperatures in the upper 12 m of the snow pack could be related to an overestimation of modelled accumulation. Pälli et al. (2002) found a relatively small accumulation rate at the top of the ice cap in comparison with the higher parts of the glacier and ascribed this difference to the effect of wind-driven snow transport. Note that the temperature at great depth is much higher than the RACMO derived yearly mean air temperature at this altitude of −12 • C, which is an indication of the significance of refreezing. In contrast to the simulated densities, the observed density profile shows a large variability between the different layers (Fig. 3). Simulating small-scale density variations would require a more detailed treatment of the evolution of the snow microstructure, melt water percolation and a higher vertical model resolution.
The aforementioned initialisation procedure requires calibrated values of the tuning parameters. Therefore, the calibration procedure in Sect. 4.4 is performed prior to the initialisation with tuned values of the calibration parameters. Since accurate calibration of γ p and C b requires initialised subsurface profiles as well, calibration experiments with various setups of γ p and C b start in 1989 (18 yr  of the actual calibration period in 2007). Note that the calibration of k, cl and b does not require initialised subsurface profiles, since the snow properties do not affect the modelled incoming shortwave and longwave fluxes.

Standard run: 1989-2010
Starting from initialised subsurface profiles, a 21-yr simulation is performed covering the period from October 1989 to October 2010. In this section, we present the temporal evolu-tion and spatial distribution of the mass and energy balance and discuss these in relation to evolving subsurface properties. Figure 4 presents contour plots of the annual mean mass balance, refreezing and runoff averaged over the full simulation period. Spatial variations in the mass balance (Fig. 4a) are to a large extent explained by a combination of altitudinal differences in air temperature and precipitation, spatial gradients in the surface albedo and local variations in the amount of refreezing. Local mass balance variations by topographic effects are only significant in the immediate vicinity of steep valley walls. The equilibrium line altitude (ELA) over the simulation period is on average 719 m a.s.l. Over the period spanned by the stake measurements (2006)(2007)(2008)(2009)(2010), the estimated observed ELA is 591 m a.s.l. This is somewhat lower than the simulated ELA over the same period of 631 m a.s.l. The spatial mean net mass balance for 1989-2010 is −0.39 m w.e. per year. In Sect. 5.3 simulated specific mass balance at the stake sites are validated against observations. As grid boundaries are selected along the ice divide, no ice is assumed to flow into the model domain. In order to compute the actual mass budget of the glacier, the negative contribution of calving should be considered in the mass budget as well. Although the glacier snout has partly retreated on land over the course of the simulation period, the mass loss by calving may still have been significant. Estimating the mass loss by calving is hindered by the absence of frontal velocity and ice thickness estimates during the active calving phase. Hence, in this study we will only discuss the surface mass balance of the glacier. Refreezing of subsurface water (Fig. 4b) is most significant in the accumulation zone, where percolating and stored water refreezes down to tens of meters below the surface. In the accumulation zone, low snow temperatures at the start of the melt season increase the potential for refreezing. Despite high melt rates in the ablation zone, refreezing is limited by the disappearance of snow during the melt season. A maximum in the amount of refreezing is found around 1000 m a.s.l. On average, refreezing contributes 0.27 m w.e. per year to the mass budget. Hence, 25 % of all melt water (1.05 m w.e. a −1 ) and rain water (0.05 m w.e. a −1 ) at the surface refreezes in the snow pack. Refreezing is equivalent to 69 % of snow accumulation during a year. The spatial pattern of runoff (Fig. 4c) is mainly controlled by the melt rate and the amount of refreezing. In the accumulation zone, a large fraction of the available water at the surface refreezes in the snow pack, thereby limiting the amount of runoff of slush water. Averaged over the glacier, the net runoff is 0.82 m w.e. a −1 , implying a total discharge of 1.58 × 10 8 m 3 a −1 .

Mass and energy budget
The surface mass balance is a product of precipitation, runoff and latent transport by sublimation and riming.   Figure 5 shows elevation profiles of the mass budget and its defining components. The mass exchange with the atmosphere by sublimation and riming is negligibly small, so the sum of precipitation and (negative) runoff determines the mass budget. High up in the accumulation zone, melt occurs while runoff goes to zero, which implies that all melt water refreezes. Precipitation dominates the mass budget at high altitudes, whereas ice melt dominates the mass budget in the ablation zone. The mean altitudinal mass balance gradient in the ablation zone is 4.1 mm w.e. a −1 m −1 , which is similar to the average found for Svalbard by Hagen et al. (2003). Contour plots of the components that comprise the energy budget are shown in Fig. 6. The sum of all the incoming and outgoing fluxes (Fig. 6a), which is equivalent to the energy involved in melting, decreases with altitude and is very small in the accumulation zone, although some melting occurs even at the highest point on the grid. The net shortwave flux (Fig. 6b) is the main source of energy, and its spatial pattern is to a large extent controlled by the surface albedo. Generally, the effect on the shortwave budget of shading near valley walls is relatively small. The net longwave flux (Fig. 6c) is the main energy sink and is strongly dependent on the temperature deficit at the surface, which increases with altitude. Regarding the turbulent fluxes (Fig. 6d), the sensible heat flux contributes significantly to the energy budget, whereas the latent heat flux is small. Due to refreezing, the subsurface heat flux (Fig. 6e) is positive and contributes significantly to the surface heat budget, especially in the accumulation zone. Due to a high surface albedo in the accumulation zone, less shortwave radiation is absorbed and surface temperatures are relatively low. This effect is partly compensated by an enhanced subsurface heat flux towards the surface. Surface heating by rainfall (Fig. 6f) is of negligible magnitude.
Time-series of the net mass balance, runoff and precipitation are shown in Fig. 7a. Over the simulation period, two positive mass balance years occurred ( 92 and 96), which are also the years with lowest runoff. The large inter-annual variability is explained by strong year-to-year variations in both summer melt and winter accumulation. The most negative net mass balance of −0.95 m w.e. a −1 is found in 1998, which is related to high summer temperatures in combination with limited snowfall during the preceding winter season. Trends in the annual and summer mean temperatures  in Fig. 7b are tested for significance by taking the ratio of regression slope and the standard error of the slope to find t-scores of 2.21 and 0.61, respectively. Consequently, only the trend in annual mean air temperatures since 1989 is significant at a 95 % confidence level. This finding is in line with observed seasonal temperature trends at Svalbard Airport since 1912, which reveal a preferred warming outside the summer months (Hanssen-Bauer et al., 2009). Correlations of −0.66 and 0.64, significant at a 99 %-confidence level, exist between net mass balance and summer mean air temperature, and between net mass balance and annual precipitation, respectively. A correlation of only −0.10 is found between the net mass balance and annual mean temperatures, which demonstrates the insignificance of temperature variations outside the melt season on the mass budget. Winter snowfall is important, owing to both its direct effect on the surface mass budget and the impact on the length of the melt season. Time-series of the annual amount of refreezing in    Figure 8a and b show time-series of the vertical temperature and density distribution for the period 1989-2010 at stake site S9, which is located in the accumulation zone (just above the ELA). Furthermore, focussed time-series for the year 2009, additionally showing the irreducible and slush water content, are presented in Fig. 8c-f. Surface melt percolates into the snow pack and subsequent refreezing heats the snow to melting point in summer. This continues until the percolating water encounters the impermeable ice and starts to form a slush layer. Heat diffuses into the underlying cold ice and after the melt season a cold wave penetrates into the firn pack from above, causing the stored slush water and irreducible water to refreeze gradually. Subsurface temperatures remain at the melting point until the water content of a layer is entirely depleted. The inter-annual variability in refreezing at S9 is quite large and determined by both surface melt water production and the cold content of the snow pack. Note that the snow depth reaches a minimum between 2004 and 2006, thereby limiting the amount of refreezing during this period. Figure 8 also shows the delayed response of ice temperatures with depth to the forcing at the snow-ice interface. An example of simulated runoff time-series at S6 in the ablation zone for the year 2010 is shown in Fig. 9. Surface melting starts in early May and until early June all melt water refreezes and runoff is absent. From then on, a slush layer is formed and water runs off gradually until the snow pack has fully melted. Clearly, the snow pack has a buffering effect on the amount of runoff; when the snow pack is all gone, all available water at the surface runs off immediately, inducing a clear diurnal cycle. A snowfall event in August shuts down ice melt and runoff for several days, which has a remarkable impact on the mass budget. A limited amount of snowfall in summer leads to a major increase in the albedo and thus effectively reduces the amount of absorbed solar radiation and thus melt.

Subsurface variables
Along the flow line marked in Fig. 1, vertical profiles of mean subsurface temperatures and densities over the final year of simulation (2009)(2010) are presented in Fig. 10. Around the ELA, a clear transition from cold ice in the ablation area to temperate snow/firn in the accumulation zone is seen. In the accumulation zone, heating by refreezing dominates advection and diffusion of lower surface temperatures at depths >5 m. In the ablation area, refreezing in the shallow snow pack during the melt season raises near-surface temperatures. Temperatures at the base of the vertical domain in the ablation zone are a measure of the mean temperature at the snow-ice interface over several decades as vertical heat diffusion is a slow process. Hence, in a changing climate this induces a buffering effect as slowly changing subsurface temperatures have a long-term impact on the subsurface heat flux at the surface. It should be noted that a zero energy flux at the lower boundary of the vertical domain might be a crude assumption in the ablation area where likely warmer ice gets advected from below. Note that inclusion of horizontal advection of ice by ice flow would lead to a horizontal shift, increasing with depth, of the vertical profiles in Fig. 10. A near-surface temperature distribution with temperate ice in the accumulation zone and cold ice near the surface in the ablation zone is typical for glaciers in Svalbard (Blatter and Hutter, 1990;Pettersson, 2004).

Comparison with observations
As discussed in Sect. 4.4, observed SW in , SW out , and LW in at the AWS between March 2009 and November 2010 are used to calibrate modelled fluxes. Scatter plots of these fluxes, presented in Fig. 11, show a good agreement between the modelled and measured values with correlation coefficients between 0.83 and 0.91. Discrepancies between observed and simulated SW in result mainly from errors in the estimated cloud cover, which also accounts for discrepancies between observed and simulated LW in . Recall that due to data-logger problems we lack observations of longwave fluxes during the melt season. The net shortwave and longwave budgets are on average overestimated by 1.6 W m −2 and 1.9 W m −2 , respectively. The only non-calibrated flux is the outgoing longwave radiation. Modelled values of LW out agree well with observations (R = 0.93), which is an indication that modelled surface temperatures are reasonably accurate. Recall that the  surface temperature is the only unknown in the energy balance equation. Validation of the surface mass balance is performed by comparing modelled values to stake observations. As mentioned, the observed surface height variations are converted into mass changes by employing the mean observed snow density in snow pits. Figure 12 reveals a good agreement between measured and simulated mass balance values for S3-11. At S1 and S2, the model overestimates mass loss, which could be related to wind-driven snow transport accumulating mass in this area. Figure 1 illustrates that, in contrast to S3-11, stakes S1 and S2 are positioned close to the southern edge of the ice grid in the vicinity of steep valley walls. Berthier et al. (2010) show that redistribution of snow by wind may lead to accumulation of snow in concave areas near the margins. Other factors that can explain the discrepancy in observed and simulated mass balance at the two sites are small scale shading and albedo effects.
Finally, we compare observed snow temperatures and densities in snow pits in spring 2008 and 2009 with simulated values. For this purpose, seven observed vertical profiles in the ablation area are used, and vertically averaged values are compared to modelled mean values. Two additional snow pit profiles were available but not used due to a lack of data points (two or less). The mean measured snow temperature of 256.2 K is only 0.5 K lower than the simulated mean. The mean measured snow density of 367 kg m −3 is 18 kg m −3 higher than simulated. Besides uncertainty in the simulated densification rate, this discrepancy may also be related to an underestimation of the fresh snow density.

Sensitivity experiments
Additional experiments are performed to test the model's sensitivity to changes in model parameters (Sect. 6.1), initial conditions (Sect. 6.2) and climate (Sect. 6.3).

Parameter sensitivity
A set of 5-yr runs over the period 2005-2010 is performed to investigate the robustness of the model output to changes in model parameters. An overview of the sensitivity of the mass balance, refreezing and runoff averaged over the glacier to perturbations of selected model parameters is given in Table 1. The simulated mass balance and runoff is most sensitive to changes in parameters that affect the net shortwave budget (α ice , α fs , k, γ p , P fs and T s/r ) and longwave radiative budget ( cl and b). The sensitivity to perturbations of the fresh snow albedo (α fs ) is particularly high, since the relative impact on the shortwave budget is large. The amount of refreezing is mainly sensitive to changes in parameters that affect the surface temperature (C b ), snow thickness (γ p ) and snow density (ρ fs ). The rate at which slush water runs off is controlled by the parameter C ro , which is a scaling factor for the runoff time-scale. Despite its impact on the runoff rate, perturbing C ro has only a limited effect on the net mass balance, refreezing and runoff. A clear nonlinear response of the mass balance to changes in the maximum accumulation rate a max is explained by a nonlinearity in the area affected by a change in a max . An upward shift of the altitude at which maximum accumulation is reached affects only a very small portion of the grid, as a relatively small amount of grid points are located above this altitude.

Initialisation sensitivity
The benefit of extensive initialisation is investigated by perturbing initial subsurface profiles. For this purpose, two simulations over the period 2004-2010 are performed with the initial water content in the snow pack set to zero in both simulations and the subsurface temperature reduced by 3 K in one of the runs. The initial water content has been set to zero to avoid immediate refreezing when lowering subsurface temperatures. Reducing the subsurface temperatures by 3 K increases the mass balance by 27 mm w.e. after 1 yr and by 60 mm w.e. between 2-6 yr after the start of the experiment. It can be concluded that lowering subsurface temperatures has only a modest but long-term effect on the mass balance, and we therefore conclude that highly-accurate simulation of the mass balance requires properly initialised subsurface conditions.

Climate sensitivity
In order to study the mass balance sensitivity to future changes in climate, we performed multiple runs with perturbed air temperature T and precipitation P over the period 2000-2010. Statistical downscaling of an ensemble of GCM future climate scenarios for Svalbard has indicated a strong seasonal variability in climate change with more pronounced warming in winter and spring (Benestad, 2008;. In one set of runs, T and P perturbations are evenly distributed over the seasons, whereas in a second set of experiments, unevenly distributed climate variations are prescribed. The degree of inhomogeneity of climate change throughout the year is based on future estimates of the seasonal mean T and P at Svalbard Airport for the period 2070-2100 . Future downscaled GCM estimates are based on the A1b-emission scenario (IPCC, 2001). Observational data from Svalbard Airport over the period 1980-2010 are used to compute mean seasonal T and P in 1995, which serve as reference values. Hence, by taking the difference of seasonal T and P between 1995 and 2085, the relative change per season with respect to the annual mean change is calculated. The related scaling factors for temperature change in winter, spring, summer and autumn are 1.48, 1.45, 0.39 and 0.68, respectively. Scaling factors for precipitation perturbations are 1.61, 1.06, 1.07 and 0.27 for winter, spring, summer and autumn, respectively. These scaling factors are used in the seasonally inhomogeneous climate sensitivity experiments (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) to scale T and P perturbations per season.
The net mass balance sensitivity to T and P perturbations with and without seasonality is presented in Table 2. Altitudinal mass balance sensitivity profiles for eight of the seasonally invariant climate change experiments are shown in Fig. 13. Clearly, the mass balance responds nonlinearly to changes in T and P . Figure 13 illustrates that the mass balance is mainly sensitive to climate variations in the ablation zone, where a lower mean surface albedo enhances the melt sensitivity. Climate perturbations lead to a shift in the ELA (Table 2) and thus the extent of the ablation zone, thereby explaining the nonlinearity. On average a 1 K temperature change can be offset by a 32 % change in precipitation of similar sign in case of seasonally uniform perturbations. On the other hand, in the case of seasonally dependent climate change, a 1 K change in temperature is balanced by a smaller change in precipitation of 17 %. Note that including seasonality has a small effect on the mass balance sensitivity to precipitation perturbations. This can be expected, since precipitation changes in winter, spring and autumn all affect the snow depth at the start of the melt season in a similar fashion. Conversely, incorporating seasonality has a major impact on the mass balance sensitivity to temperature perturbations. Since the mass balance is mainly a product of melt in the summer months and year-round precipitation, a relatively small increase of summer temperatures has a pronounced The Cryosphere, 6, 641-659, 2012 www.the-cryosphere.net/6/641/2012/ effect on the mass balance sensitivity. This stresses the importance of considering climate change on a seasonal scale, especially in regions like the Arctic where temperature increase is expected to be less pronounced during the melt season. Another remarkable feature in Table 2 is that the amount of refreezing is not very sensitive to changes in both temperature and precipitation. An increase in the amount of precipitation leads to an increase in the amount of refreezing, since a thicker snow pack at the start of the melt season potentially has a larger cold content. Regarding the sensitivity of refreezing to temperature changes, two effects seem to counteract each other: (1) the accumulation zone extent decreases with increasing temperature, which reduces the potential for refreezing, and (2) the surface melt water production increases with increasing temperature, thereby enhancing the potential for refreezing. Seasonally dependent sensitivities in Table 2 indicate that the amount of refreezing does not change significantly in a future climate for the fixed geometrical setup of the sensitivity experiments. In order to put the mass balance period 1990-2010 into a longer term perspective, we first used output of additional runs to construct Seasonal Sensitivity Characteristics (Oerlemans and Reichert, 2000). In these model experiments, running over the period 2000-2010, T or P are perturbed for one season, while unperturbed values are used for the other Table 2. Overview of the mass balance sensitivity (δMB), refreezing sensitivity (δRE) and ELA sensitivity (δELA) to changes in air temperature T and precipitation P . Sensitivities are given in m w.e. a −1 for the mass balance and refreezing sensitivity and m a.s.l. for the ELA sensitivity. The non-seasonal δELA in the run with T + 4 K could not be resolved since the ELA in this scenario is above the highest point on the grid. Mass balance, refreezing and ELA values in the unperturbed climate run are −0.448 m w.e. a −1 , 0.259 m w.e. a −1 and 724 m a.s.l., respectively.

Run
No  seasons. In this way the sensitivity of the annual mass balance to seasonal climate variations can be estimated, and the sensitivities will not depend on changes during the other seasons (affecting snow depth and subsurface conditions). Figure 14 shows the mass balance sensitivity to various seasonal T and P changes. The mass balance sensitivity to T variations is pronounced in summer, whereas in winter the impact is negligible. Seasonal P perturbations have a yearround impact on the annual mass balance. For the period 1912-1989, we derived time-series of seasonal mean T and P values from monthly composite time-series at Svalbard Airport . Mean seasonal temperatures and precipitation totals were computed for the period 1912-1989 and 1990-2010 and indicate an increase in annual mean temperature and annual precipitation of 1.8 • C and 11 mm (6.2 %), respectively. When looking at the seasonal pattern in temperature change, in summer the increase is much less pronounced (1.0 • C) than in winter (2.3 • C), which is in line with the aforementioned future projection for this region. Nevertheless, from Fig. 14 it can be concluded that the modest increase in summer temperatures had a larger impact on the mass balance than strongly rising winter temperatures. Based on the trends in temperature and precipitation and sensitivities presented in Fig. 14, we conclude that the mass balance was more positive during the period 1912-1989 than during the simulation period . Note that using the annual mean change in T and P and seasonally invariant sensitivies (Table 2) to estimate the trend in the mass balance would lead to an overestimation. We purposely do not present absolute values or time-series of the mass balance over the period 1912-1989 or in the future, since ice dynamical changes affect the ice geometry (glacier area and surface height) and lead to an increasingly large uncertainty in mass balance estimates over longer periods of time. Factors not taken into account include the mass balance altitude feedback, changes in glacier length and area, changes in grid orientation (affecting absorbance of solar radiation) and changes in the subglacial hydrological system (affecting the amount of basal lubrication). Finally, the mass balance effect of long-term trends (>10 yr) in subsurface conditions is not captured in the relatively short period of the sensitivity experiments. Estimating the magnitude of the mass balance effect of these factors requires coupling of the mass balance model to an ice flow model, which is left for future work.

Conclusions and discussion
In an attempt to simulate the surface mass budget of Nordenskiöldbreen, including refreezing and runoff, a distributed surface energy balance model is coupled to a multi-layer snow model. Gridded climate input is generated from regional climate model (RACMO) output in combination with meteorological time-series from Svalbard Airport. An optimal parameter setup is found after calibration with observations on the glacier and at Svalbard Airport. Prior to the actual simulation over the period 1989-2010, extensive initialisation of subsurface profiles is performed to find initial - subsurface conditions that evolve to observed borehole profiles in 1997.
Over the period 1989-2010, the simulated net surface mass balance is −0.39 m w.e. a −1 with extrema of +0.17 and −0.95 m w.e. a −1 for mass balance years 91-92 and 97-98, respectively. Year-to-year variations are explained by fluctuations in summer air temperatures and yearround snowfall. Refreezing of subsurface water amounts to 0.27 m w.e. a −1 , which is equivalent to 69 % of the annual snow accumulation. The mean simulated ELA is 719 m a.s.l. The net runoff per unit area is the sum of discharge of slush water and surface runoff and is equal to 0.82 m w.e. a −1 . Refreezing delays and reduces runoff especially in the accumulation zone, where refreezing is most pronounced. The simulated thermodynamic structure is characterised by a temperate accumulation zone (at depths >∼5 m) and cold ice in the ablation zone. Large amounts of refreezing in the accumulation zone are a consequence of low annual mean temperatures (<−10 • C in the accumulation zone), providing a large cold content in the snow/firn pack at the start of the melt season, in combination with significant melt during the summer season. Conversely, in the ablation area annual snow accumulation mainly controls the amount of refreezing.
Modelled energy fluxes are validated against AWS measurements on the glacier, and a generally good agreement is found. A comparison of mass balance estimates from stake observations to simulated values shows a good agreement, except for the two lowest stake locations. We hypothesise that this is related to the effect of wind driven snow transport not considered by the model. Snow pit profiles in the ablation zone indicate that simulated snow temperatures correlate well with observations, whereas the density is somewhat underestimated, which is shown to have a very limited impact on the mass balance.
Parameter sensitivity experiments show a high sensitivity of the mass balance and runoff to parameters affecting the shortwave and longwave radiative budget. The amount of refreezing is sensitive to changes in parameters that affect the surface temperature, snow thickness and snow density. Experiments with perturbed initial subsurface conditions show that the simulated mass balance is modestly affected by changes in initial subsurface temperatures.
Multiple climate perturbation experiments are performed and show a non-linear response of the mass balance, refreezing and runoff to changes in temperature and precipitation. Including seasonal variations in climate change, based on a future climate scenario for Svalbard Airport, reduces the mass balance sensitivity, mainly as a result of less pronounced summer warming. Seasonal sensitivity characteristics confirm that the mass balance is particularly sensitive to temperature changes in summer, whereas the seasonal variability in mass balance sensitivity to precipitation changes is relatively small. Temperature and precipitation time-series from Svalbard Airport since 1912 reveal colder and somewhat drier conditions over the period 1912-1989 with respect to the period 1990-2010, which makes it plausible that the period preceding the simulation period had a more positive mass balance. In line with future projections for this region, these time-series also reveal least pronounced warming during the melt season. With the fixed geometrical setup of the climate perturbation experiments, due to compensating effects, refreezing is hardly affected in a future climate.
Many factors lead to uncertainties in the model output, of which probably the largest source of error is in the model input. One of the main uncertainties in the model forcing is related to spatial and temporal variability in snow accumulation. In the model, the accumulation rate increases linearly with height (at a rate γ p ) and maximizes at a max . In reality, wind patterns will have a major impact on the spatial distribution of the accumulation rate. As discussed by Berthier et al. (2010), a redistribution of snow by wind often leads to enhanced thinning rates along the centerline and additional accumulation in concave areas near the margins. Additionally, temporal variability in γ p will also have a considerable impact on the simulated mass balance. Accumulation rates along a horizontal profile in the accumulation zone of Nordenskiöldbreen, presented by Pälli et al. (2002), confirm that the accumulation rate varies significantly in space. Errors in the cloud cover estimates, derived from observations at Svalbard Airport, have a significant impact on the modelled shortwave and longwave budget. A comparison of observed air temperatures at the AWS site and RACMO derived temperatures shows a high correlation, and systematic errors for this specific site are found to be small. In the absence of sufficient meteorological time-series near the glacier, output of RACMO therefore provides valuable data to force the model. Uncertainties in the simulated mass balance distribution are additionally related to unmodelled effects, like wind-driven snow drift and dust deposition. Regarding the subsurface profiles uncertainties arise from the zero energy flux assumption at the base of the vertical domain, the limited amount of subsurface layers and disregarding horizontal advection of heat, mass and water.
Finally, it is worth mentioning that output fields of the subsurface evolution may serve as surface boundary input for an ice dynamical model. Realistically simulating the geometric evolution of a glacier requires dynamic coupling of an ice flow model to a surface model that provides details on the water input, surface temperature and the mass balance. Basal sliding velocities at the ice-bedrock interface are known to depend strongly on the rate of water input from the surface (Zwally et al., 2002;Van de Wal et al., 2008;Schoof, 2010). Refreezing in the snow and firn pack at the start of the melt season effectively delays water transport to the base and hence affects the timing of speed-up events, as observed during the onset of the melt season on Nordenskiöldbreen (Den Ouden et al., 2010). This demonstrates the importance of incorporating a surface mass balance model that accounts for refreezing in order to model ice velocities at seasonal time-scales and smaller.