Modelling the long-term mass balance and firn evolution of glaciers around Kongsfjorden, Svalbard

Abstract

We analyse the long-term (1961–2012) distributed surface mass balance and firn evolution of the Kongsvegen and Holtedahlfonna glacier systems in northwestern Svalbard. We couple a surface energy-balance model to a firn model, with forcing provided from regional climate model output. In situ observational data are used to calibrate model parameters and validate the output. The simulated area-averaged surface mass balance for 1961–2012 is slightly positive (0.08 mw.e.a⁻¹), which only fractionally compensates for mass loss by calving. Refreezing of percolating water in spring/summer (0.13 m w.e. a⁻¹) and stored water in fall/winter (0.18 m w.e. a⁻¹) provides a buffer for runoff. Internal accumulation, i.e. refreezing below the previous year’s summer surface in the accumulation zone, peaks up to 0.22 m w.e. a⁻¹, and is unaccounted for by stake observations. Superimposed ice formation in the lower accumulation zone ranges as high as 0.25 m w.e. a⁻¹. A comparison of the periods 1961–99 and 2000–12 reveals 21% higher annual melt rates since 2000 and a 31% increase in runoff, which can only in part be ascribed to recent warmer and drier conditions. In response to firn line retreat, both albedo lowering (snow/ice–albedo feedback) and lower refreezing rates (refreezing feedback) further amplified runoff.

1. Introduction

Despite the relatively small ice volume of Svalbard glaciers, compared to the ice sheets of Greenland and Antarctica, their high sensitivity to rapidly changing Arctic climate has led to increasing interest in studying their evolution. The coupling between sea-ice retreat around the archipelago, the surface reflectivity and air temperature is thought to have caused amplified warming (‘Arctic amplification’) with respect to the global mean (Serreze and Francis, 2006; Serreze and others, 2009). Temperature in Longyearbyen, Svalbard, has risen by 0.27 K decade⁻¹ since 1912, with more pronounced warming in winter/spring (0.38 K decade⁻¹) than in summer/autumn (0.16 K decade⁻¹). These trends are expected to continue in the 21st century, yielding projected annual mean temperature changes between 1961–90 and 2071–2100 ranging from 3 K in the southwest to 8 K in the northeast of Svalbard (Førland and others, 2011).

The interaction between a glacier surface, the atmosphere above and the underlying firn determines the surface mass balance, defined as the sum of mass gain and loss per unit area. While stake observations in combination with density measurements provide in situ point estimates of the seasonal (summer and winter) balance, mass-balance models are useful tools for simulating spatio-temporal variability of the mass balance and its contributing components in greater detail and over longer periods of time. Detailed modelling of the mass balance requires solving the surface energy balance to obtain the surface temperature and melt production, and then coupling the surface model to a subsurface routine to account for the impact of firn water storage and refreezing on the mass- and energy budget.

Energy consumed in melting the surface may be released again in the snowpack after refreezing of percolating or stored water. Refreezing not only acts to heat the snow, it also acts as a buffer for mass loss since part of the melt produced at the surface is retained in the snowpack, thereby increasing the snow density. Its significance for the mass budget of polythermal Arctic glaciers is known to be substantial (Wright, 2005; Reijmer and Hock, 2008; Van Pelt and others, 2012). Refreezing has a seasonal cycle and is generally most pronounced at the start of the melt season, when melt first enters the cold snowpack (Schneider and Jansson, 2004). In the course of the melt season, refreezing occurs preferentially near the surface during melt–freeze cycles (Colbeck, 1986) and may be limited by the absence of a snowpack in the ablation area. After the melt season, refreezing of stored water continues in the accumulation area as the winter cold wave penetrates the snow and firn (Pfeffer and others, 1991; Schneider and Jansson, 2004).

When refreezing in the accumulation zone occurs below the previous year’s summer surface, this is generally referred to as internal accumulation. As stake mass-balance measurements only account for mass changes above the previous year’s summer surface, internal accumulation may lead to a substantial overestimation of mass loss from stake observations in the accumulation zone (Trabant and Mayo, 1985; Schneider and Jansson, 2004; Reijmer and Hock, 2008). At the subsurface interface from snow or firn to ice, water tends to collect and it may refreeze when cooling occurs from above or below to form an ice layer. When this refrozen ice layer is exposed at the surface it is generally referred to as superimposed ice (SI). At the end of the melt season, SI is found just above the equilibrium line in the superimposed zone and may complicate visual detection of the equilibrium-line altitude (ELA). As SI formation may constitute a substantial component of the mass budget (e.g. Wadham and Nuttall, 2002; Bøggild and others, 2005; Wright and others, 2007; Brandt and others, 2008), it is important to account for when converting stake height measurements into mass-balance estimates.

Predicting mass-balance evolution in a changing climate is mostly a matter of understanding the role of feedbacks between atmospheric, surface and firn processes. A wellknown example is the height/mass-balance feedback, in which enhanced melt lowers the surface height, thereby increasing air temperatures, which further accelerates melt (e.g. Oerlemans, 2001). Another significant feedback is the snow–albedo feedback, which arises when enhanced melt leads to a retreat of the firn line, increased bare ice exposure and a drop in the mean albedo, which further enhances surface melt (e.g. Hall, 2004; Van den Broeke and others, 2008). The mass-buffering effect of refreezing tends to slow down these feedback processes, but its dampening effect is likely to be reduced in a warming climate, implying accelerated glacier mass loss (Wright and others, 2007; Van Angelen and others, 2013).

In this study, we couple a surface energy-balance and snow–firn model to simulate the long-term distributed mass balance and firn evolution of glaciers around Kongsfjorden in northwestern Svalbard. The surface energy-balance model is forced with downscaled meteorological output from the HIRLAM regional climate model, covering the period 1961–2012 (Undén and others, 2002; Reistad and others, 2011). The subsurface model simulates the vertical evolution of density, temperature and water content, while accounting for percolation, refreezing, storage and runoff of water. A wealth of observational in situ data from automatic weather stations (AWSs), stake measurements and shallow firn cores are used to calibrate parameters and validate the model output. With this work we aim to: (1) estimate and analyse the long-term (1961–2012) evolution of the distributed surface mass balance; (2) quantify refreezing, internal accumulation, and SI formation; and (3) detect multi-decadal trends in surface and subsurface conditions by comparing the periods 1961–99 and 2000–12.

2. Study Area and Previous Work

The focus of this study is on five glaciers around Kongsfjorden, shown in Figure 1: Sidevegen, Kongsvegen, Infantfonna, Fatumbreen and Holtedahlfonna/Kronebreen.

Fig. 1. Surface-height map of the study area including outlines of the glaciers Sidevegen, Kongsvegen, Infantfonna, Fatumbreen and Holtedahlfonna/Kronebreen. The inset map shows the location of the glaciers in western Svalbard. Black dots mark the position of mass-balance stake measurements; orange circles mark the positions of shallow core sites; blue circles mark the locations of AWSs. Red crosses mark the midpoints of the 11 km × 11 km grid of the regional climate model, used to generate the downscaled 100 m × 100 m climate forcing for the coupled model (Section 3.3).

The five glaciers have individual accumulation basins, but all merge downstream and enter the fjord with a single front. We distinguish between two glacier systems: the Holtedahlfonna system (HDF), consisting of Holtedahlfonna/Kronebreen, Infantfonna and Fatumbreen, and the Kongsvegen system (KNG), comprising Kongsvegen and Sidevegen. Based on a digital elevation model (DEM) from 2009 (NPI, unpublished data), the Kongsvegen and Holtedahlfonna systems span a total area of 173 and 385 km² and cover an altitudinal range of 0–1059 and 0–1441 m a.s.l. respectively. In terms of ice dynamics the two glacier systems are very different. Kongsvegen is a surge-type glacier, which last surged in 1948 (Melvold and Hagen, 1998), and which is currently in its quiescent phase, with low flow velocities (2–8 m a⁻¹) and surface steepening (Hagen and others, 1999; Nuth and others, 2012). Conversely, the lower reach of HDF, Kronebreen, is a steady fast-flow tidewater glacier with ice velocities in the range 300–800 m a⁻¹ (Kääb and others, 2005). Between 1966 and the late 1980s, the front of Kronebreen retreated ∼3 km. It remained at more or less the same position from 1990 until 2011, when rapid frontal retreat began, with the grounded calving front retreating ∼1 km between 2011 and 2013.

Several previous studies have examined the energy- and mass balance of glaciers around Kongsfjorden. The distributed surface mass balance of KNG and HDF was investigated by Nuth and others (2012), who applied a simple degree-day model to obtain surface melt rates, used as part of a mass continuity method to estimate glacier mass loss related to calving. Nuth and others (2012) found a weakly negative net surface mass balance for KNG and HDF of −0.04 and −0.02 m w.e. a⁻¹ between 1969 and 2007, with accelerated mass loss since 1990. In earlier studies, Lefauconnier and others (1994) used radioactive layering to estimate a slightly negative net mass balance (−0.11 m w.e. a⁻¹) for Kongsvegen over the period 1967–90, while Rasmussen and Kohler (2007) used a simple statistical approach to reconstruct the mass balance of Kongsvegen back to 1948. Most recently, Karner and others (2013) investigated the energy- and mass balance for 2001–10 at a single site on Kongsvegen, using a coupled surface–snow model, driven by meteorological data from an AWS. They found a rough balance between net solar and thermal radiation and detected the formation of SI during positive mass-balance years. SI formation had earlier been detected from satellite imagery on Kongsvegen by König and others (2002), and Obleitner and Lehning (2004) successfully reproduced observed SI with a physically based snow model. Additionally, a radar survey in the lower accumulation zone of Kongsvegen by Brandt and others (2008) revealed substantial amounts of SI (up to 0.43 m w.e. a⁻¹), formed since the 1960s.

3. Model and Set-Up

We apply a coupled surface-energy-balance–snow model to simulate mass- and energy exchange between the atmosphere, surface and underlying snow, firn and/or ice. The coupled model has previously been used by Van Pelt and others (2012) to analyse the mass balance and climate sensitivity of Nordenskiöldbreen, a tidewater glacier in central Svalbard. More recently, the coupled model was applied in inverse approaches to Nordenskiöldbreen to reconstruct basal topography (Van Pelt and others, 2013) and to determine spatio-temporal accumulation variability from ground-penetrating radar data (Van Pelt and others, 2014).

3.1. Model description

Here, only the main qualitative aspects of the model are outlined. Van Pelt and others (2012) have published a more complete overview of the coupled model and its implementation.

The surface energy-balance model computes all the energy fluxes at the surface interface and determines the surface temperature for which the net energy flux is zero, i.e. for which energy is conserved. Surface temperatures are held at the melting point and excess energy is used for melting. The energy fluxes contributing to the energy budget comprise incoming and reflected shortwave radiation, incoming and emitted longwave radiation, sensible and latent turbulent heat exchange, heat supplied by rain and the subsurface heat flux.

Energy- and mass exchange at the surface provides upper boundary conditions for the subsurface model, tracking the subsurface evolution of temperature, density and water content. A ‘tipping-bucket method’ is used to account for vertical water transport, storage and runoff (Greuell and Konzelmann, 1994; Reijmer and Hock, 2008; Van Pelt and others, 2012). Water may be held in the snow against capillary forces as ‘irreducible water’ or it may be stored in remaining pore spaces as ‘slush water’ on top of an impermeable ice layer. Refreezing of percolating or stored water in firn or snow occurs when subsurface temperatures are below melting point. Refreezing not only adds solid mass to a snow layer, increasing its density, it also releases heat, raising the layer’s temperature. As the energy consumed in melting the surface may be released again below the surface, meltwater percolation and refreezing can be regarded as a form of latent heat transport. Runoff is defined as the amount of vertical water exchange at the firn–ice transition and occurs instantaneously in the case of bare-ice exposure and according to a slope-dependent timescale for slush water in the firn/snowpack (Reijmer and Hock, 2008). No horizontal exchange of slush water occurs between gridcells, and local runoff follows the principles of a linear reservoir model. The subsurface model solves the thermodynamic equation, accounting for heat conduction/ diffusion and refreezing, and it solves the densification equation, considering refreezing and gravitational packing of the snow. Vertical advection of mass and energy is accounted for by implementing a moving grid, which may add/remove layers at the surface and base of the vertical domain, in response to surface height changes (Section 3.2).

The surface mass balance is defined as the accumulated mass exchange per unit area and is the sum of mass gain by precipitation and condensation or riming, and mass loss by runoff and evaporation or sublimation. Refreezing acts as a buffer for mass loss as it may reduce and delay runoff where a firn/snowpack is present. The surface mass balance is therefore effectively a product of surface and subsurface processes. For calibration of model parameters against stake data (Section 3.4), it is useful to also compute the ‘stake mass balance’. In contrast to the actual mass balance defined above, the simulated stake mass balance only accounts for mass changes above the previous year’s summer surface and is therefore equivalent to what is observed with mass-balance stakes. The difference between the actual and stake mass balance is equal to the amount of water that percolates past the previous year’s summer surface and is retained deeper in the firn pack in the form of internal accumulation or as stored water.

3.2. Model grid

Two DEMs are used to construct a time-dependent elevation grid for the model experiments. A first DEM, with a 10 m spatial resolution, was constructed by the NPI from aerial photogrammetric imagery, taken in 1966, in combination with contour data to fill in remaining gaps (Altena, 2008). A second DEM, with a 5m horizontal resolution, was constructed from aerial photogrammetric imagery taken in 2009 (Fig. 1). Averaging was applied to both DEMs to obtain 100 m × 100 m grids with a total of 157 491 gridpoints. By differencing the DEMs and dividing by the 43 year period between the acquisition dates, we determined a distributed pattern of the linear rate of change of surface elevation in time between 1966 and 2009. This rate-of-change pattern is then applied to the full simulation period (1961–2012), to obtain an evolving surface geometry with a gradual transition from the 1966 to the 2009 surface height pattern. The above provides an efficient alternative to the more sophisticated approach of explicitly modelling ice dynamics coupled to surface conditions. Outlines for the KNG and HDF glacier systems were taken from the GLIMS (Global Land Ice Measurements from Space) database (König and others, 2013; Nuth and others, 2013), reflecting the 2007 extent of the glaciers. We chose to apply fixed outlines for the entire simulation period to allow for comparison of spatially averaged fields and trends over the same area. Hence, spatial patterns are computed for a reference surface with fixed extent and linearly varying altitude in time.

The subsurface model uses a moving vertical grid with a total of 150 layers. The uppermost model layer has a variable depth, and model layers shift downwards when the upper layer thickness exceeds 0.1 m, or shift upwards if the upper layer vanishes. Layer thicknesses are increased at depth from the surface by merging at the 25th, 50th and 75th layers. Gravitational densification reduces layer depths as snow or firn layers age. At the base of the vertical domain a Neumann boundary condition is applied, in which we assume that the temperature and density gradient between the two lowermost layers also applies at the base.

3.3. Climate forcing

As a meteorological forcing, the coupled model requires meteorological input of air temperature, humidity, cloud cover, air pressure and precipitation on a distributed grid. As a main source, we use regionally downscaled European Centre for Medium-Range Weather Forecasts reanalysis data (Uppala and others, 2005), obtained using the hydrostatic numerical weather prediction model High Resolution Limited Area Model (HIRLAM; Undén and others, 2002; see www.hirlam.org for more documentation), which yields 3 hourly meteorological fields at an 11 km horizontal resolution (Reistad and others, 2011). Midpoint locations of the regional climate model (RCM) data are marked in Figure 1.

Further downscaling is done to project the RCM output onto the 100 m × 100 m coupled model grid. Different downscaling approaches are applied per climate variable depending on the significance of topography on the distributions and the availability of observational data for calibration. Three-hourly specific humidity and cloud cover output at the RCM gridpoints are linearly interpolated onto the 100 m × 100 m grid without further corrections. In order to account for subgrid orographic effects, RCM air temperature and pressure fields are first projected to sea level, using lapse rates determined by fitting a linear (temperature) and exponential (pressure) function to the RCM data points and their known elevations. In a next step, sea-level temperature and pressure fields are spatially interpolated onto the 100 m × 100 m grid, and then projected to the DEM elevations using the previously determined lapse rate functions. To generate precipitation input, spatial mean values of all the RCM gridpoints in the domain are computed at every 3 hour time slice. The resulting RCM precipitation time series are projected onto the model grid using elevation-dependent functions, which are calibrated against observational data from centre-line winter stake mass-balance measurements (further discussed in Section 3.4). Any possible systematic biases in the temperature, pressure, humidity and cloud-cover fields are not directly corrected for, but are indirectly compensated for during calibration of surface fluxes and stake balances (Section 3.4).

3.4. Parameter estimation

3.4.1. In situ measurements

In order to calibrate model parameters, we use in situ measurements from mass-balance stakes, snow density pits, snow depth probing and AWSs.

Exposed stake-height changes along the centre lines of Kongsvegen and Holtedahlfonna have been measured biannually in April/May and September (Hagen and others, 1999, 2005). On Kongsvegen (K1–9), the stake mass-balance record starts in 1987, while seasonal mass-balance observations on Holtedahlfonna (H1–10) started in 2003. On the rough crevassed surface of the lower part of Holtedahlfonna (Kronebreen; KR1–3), mass-balance measurements have been performed only since 2008 due to limited accessibility of the terrain. Stake height readings are combined with snowpack density measurements from snow pits at selected sites every spring to yield winter mass-balance estimates, and summer balance is calculated by the change in vertical water equivalent between spring and autumn (Hagen and others, 1999; Nuth and others, 2012).

Meteorological data from two AWSs at K6 on Kongsvegen (since May 2007) and H4.5 on Holtedahlfonna (since May 2009) include flux measurements of incoming and outgoing shortwave and longwave radiation, measured using a Kipp & Zonen CNR1 radiometer. The two AWSs are located at altitudes of 534 (K6) and 684 m a.s.l. (H4.5) (Fig. 1). Processing of the data comprised removal of outliers related to logger failure and sensor riming, and data are averaged to obtain 3 hourly values. No corrections for possible misalignment of the sensors are applied.

3.4.2. Calibrating parameters

Following a similar procedure to that of Van Pelt and others (2012), we calibrate parameters determining the incoming shortwave radiation, incoming longwave radiation, snow/ firn albedo, precipitation distribution and summer melt. We calibrate the parameters in sequence, with the order of calibration chosen carefully to reduce ambiguity in the constrained parameter values. The order is determined by considering the dependence of quantities to be calibrated on the calibration of the other parameters. Hence, we start by calibrating the independent quantities such as the incoming radiative fluxes, which are unaffected by the calibration of the other calibration constants afterwards. The final parameters to be calibrated determine summer melt, which depends most strongly on the selected values of the other calibration constants.

Model calibration experiments run over the observation period for the gridpoints of interest, after initializing the model for a short period of 5 years. Longer initialization is irrelevant for calibrating parameters only affecting near-surface processes, but is required prior to the full model run, as discussed in Section 3.5. Table 1 gives an overview of the calibration parameters and the related model quantities that are matched to the observations.

Table 1. Overview of calibration parameters. Columns list, respectively, the calibration parameters, the symbol used, calibration value and unit, the matched quantity used to optimize the calibration parameter, and the observation sites used for calibration (Fig. 1)

The first calibration parameters (k, ∊ _cl and b) are set by matching modelled and observed mean incoming shortwave radiation, incoming longwave radiation under clear-sky conditions and incoming longwave radiation under cloudy conditions, respectively (Van Pelt and others, 2012). Scatter plots in Figure 2a and b show simulated vs observed incoming shortwave and longwave fluxes at K6 and H4.5. No substantial offsets are apparent between the two locations, so we decided to select fixed values for k (0.974), ∊ _cl (0.455) and b (0.960) for the entire grid (Table 1). Values for k, ∊ _cl and b are found to agree well with values reported by Klok and Oerlemans (2002) and Van Pelt and others (2012).

Fig. 2. Comparison of simulated and observed 3 hourly incoming shortwave radiation (a) and incoming longwave radiation (b), observed with AWSs at H4.5 and K6. Black lines represent the reference 1 : 1 line.

Calibrated values for the fresh snow albedo α _fs (0.80) and firn albedo α _firn (0.52) are selected by matching the mean maximum (after snowfall) and minimum (old snow) modelled surface albedo at K6 and H4.5 to observed albedo estimates, obtained by taking the ratio of reflected and incoming solar radiation.

The RCM output provides 3 hourly time series of spatial mean precipitation (Section 3.3). To generate distributed spatial patterns of precipitation for every model time step, precipitation–altitude functions are constructed and coefficients are calibrated to reduce the misfit between observed winter balance measurements and simulated values at the stake sites on KNG and HDF. In congruence with the observed winter balance, the modelled winter balance in spring represents the accumulated mass on top of the previous year’s summer surface. As in Van Pelt and others (2012), precipitation can either fall as rain or snow and a gradual transition occurs around a temperature threshold of 1.5°C. For the precipitation–altitude relation, a best match between modelled and observed winter balances is achieved when assuming precipitation P(z, t) increases linearly with elevation:

(1)

where P _RCM(t) is the RCM mean precipitation time series, z _RCM is the mean altitude of the RCM gridpoints (602ma.s.l.), z is the grid altitude, and C ₁ and C ₂ are coefficients, used to calibrate the mean precipitation and altitudinal gradient respectively. C ₁ and C ₂ are calibrated to optimize agreement between mean altitudinal profiles of stake winter mass balance, averaged over the full observation period. Like Nuth and others (2012), we found a steeper precipitation gradient along the centre line of Kongsvegen, compared to Holtedahlfonna; this led us to split the domain and determine independent values for C ₁ and C ₂ for the KNG and HDF system (Table 1). Precipitation on HDF is found to reach a maximum in the upper accumulation area, above which it is approximately constant. In the model we hold precipitation on HDF at a constant level above 1000 m a.s.l., but no upper limit has been set on KNG. The two precipitation–altitude functions are found to intersect at 160 m a.s.l., which is about the altitude at which the HDF and KNG systems merge. As the merged part of KNG and HDF is likely to experience a similar precipitation forcing, the better-constrained KNG forcing is applied to the HDF system in its lower part. After calibration of the precipitation forcing, modelled and observed winter stake balances agree well (R = 0.91) for all stakes (Fig. 3a and b). Time series of the mean winter balance at KNG and HDF in Figure 3c and d demonstrate that the model is well able to reproduce year-to-year variations, thereby giving confidence in the interannual variability generated by the RCM.

Fig. 3. Comparison of simulated and observed summer, winter and net stake mass balances. Scatter plots in (a) and (b) compare modelled and observed stake balances for all stake sites on HDF (a) and KNG (b) during the full period of observations. Time series in (c) and (d) compare modelled and observed centre-line averaged stake balances for HDF (c) and KNG (d).

In the last step of the calibration procedure, the albedo of ice α _ice (0.39) and the turbulent exchange coefficient C_b (0.0025) are collectively tuned to reduce the misfit between the long-term mean modelled and observed summer stake mass balances (Table 1). The two parameters could uniquely be determined to minimize the summer balance misfit as α _ice becomes increasingly important for the summer balance towards lower altitudes, whereas C_b has a glacier-wide impact on the summer balance. A high correlation between modelled and observed summer balances is found (R = 0.95), as illustrated in Figure 3a and b. Furthermore, interannual variability in the summer mass balance is well captured by the model (Fig. 3c and d).

After calibration, simulated net stake mass balances are found to correlate very well with the observations (R = 0.96; Fig. 3). Remaining discrepancies (RMSE of 0.24 m w.e.) can be ascribed to uncertainties in the observations, model physics and climate forcing. As for the winter and summer stake balance, time series of the stake-averaged net mass balance are in good agreement, with the model capturing most of the observed interannual variations (Fig. 3c and d).

3.5. Initialization

Prior to the full model run (1961–2012), an initialization run is done to generate initialized subsurface conditions (temperature, density and water content) at the start of the simulation. In the absence of a long record of climate input prior to 1961, initialization consisted of running the model twice over the period 1961–90 (starting the second iteration from the final state of the first iteration). At the start of the first iteration, the vertical grid consisted of ice with a constant temperature of 270 K everywhere. After two loops of initialization, subsurface conditions during the full run in 1990 turned out to be nearly identical to subsurface conditions at the end of the initialization, thereby indicating sufficient spin-up of the model. As there is uncertainty in the representativeness of the climate forcing during initialization, sensitivity tests with perturbed climate during initialization were performed, and are discussed in Section 4.2.

4. Results and Discussion

After preparing the climate input (Section 3.3), setting model calibration constants (Section 3.4) and initializing the model (Section 3.5), a 51 year simulation covering the period 1961–2012 is performed.

4.1. Mass balance, refreezing and internal accumulation

4.1.1. Spatial variability

Figure 4 shows the long-term (1961–2012) mean distributed fields of the surface mass balance, subsurface refreezing and internal accumulation. Averaged over the entire grid, the surface mass balance is found to be slightly positive (0.08 m w.e. a⁻¹). For KNG and HDF we find gridded mean surface mass-balance estimates of 0.01 and 0.13 m w.e. a⁻¹, respectively. From mass continuity, Nuth and others (2012) estimated calving on KNG to have been insignificant over the simulation period, so the surface mass balance can be regarded as a direct measure for glacier mass changes. In response to the steep mass-balance–altitude gradient and the weak mass turnover (Melvold and Hagen, 1998), the surface topography of Kongsvegen is steepening and the glacier is likely to build up for a new surge. For HDF, the positive mean mass balance (0.13 m w.e. a⁻¹) is not nearly sufficient to balance mass losses by calving of −0.37 to −0.52 m w.e. a⁻¹, estimated by Nuth and others (2012), implying substantial glacier thinning over the simulation period. The ELA is lower on KNG (500 m a.s.l.) than on HDF (610 m a.s.l.), which can be ascribed to the steeper precipitation gradient on KNG, implying higher precipitation rates up-glacier. The surface mass balance is found to be dominated by mass gain due to precipitation (0.87 m w.e. a⁻¹) and mass loss due to runoff (0.79 m w.e. a⁻¹); net latent exchange of mass is found to be a small term (0.01 m w.e. a⁻¹). Topographic effects (e.g. shading by surrounding topography or the surface slope orientation) have some impact on the mass balance, particularly in the vicinity of valley walls, although the overall influence is small.

Fig. 4. Long-term mean spatially distributed patterns of the mass balance (a), refreezing (b) and internal accumulation (c), averaged over the simulation period 1961–2012.

Refreezing of percolating and stored melt or rain in snow and firn amounts to 0.30 m w.e. a⁻¹ on average (Fig. 4b), and is similar to the mean value of 0.27 m w.e. a⁻¹ found on Nordenskiöldbreen, central Svalbard (Van Pelt and others, 2012). The amount of refreezing is equivalent to 35% of the net annual precipitation and is therefore a significant component of the surface mass budget. Figure 4b illustrates that refreezing increases with altitude and that it is most pronounced in the accumulation area.

A substantial fraction of refreezing in the accumulation zone occurs below the summer surface of the previous year, yielding an area-averaged mean internal accumulation rate of 0.11 m w.e. a⁻¹, which is equivalent to 37% of total refreezing (Fig. 4c). On average, internal accumulation comprises 13% of annual precipitation. Locally, internal accumulation amounts to up to 0.22 m w.e. a⁻¹. These high internal accumulation rates suggest a substantial underestimation of the mass balance when measured with stakes in the accumulation zone. Recall that this observational bias does not induce a bias in the modelled mass balance after calibration, as we used the simulated ‘stake mass balance’, ignoring the effect of internal accumulation, rather than the actual mass balance, for parameter optimization, as discussed in Section 3.1. To our knowledge, we present the first estimates of internal accumulation on glaciers in Svalbard. Previous studies indicate substantial amounts of internal accumulation on Storglaciären, Sweden, ranging between 3% and 5% (Schneider and Jansson, 2004) and 20% (Reijmer and Hock, 2008) of annual accumulation, while, for five Alaskan glaciers, Trabant and Mayo (1985) quantified internal accumulation at 7–64% of annual accumulation.

Refreezing can be split into three components: (1) refreezing of percolating water, (2) refreezing of stored irreducible water and (3) refreezing of stored slush water. The long-term mean spatial patterns of these three components are shown in Figure 5a–c. Refreezing of percolating water (Fig. 5a) is most pronounced in spring, when the first melt enters the cold snow or firn and refreezes. Additionally, rainfall during winter or spring can be a substantial component. Refreezing of percolating water reaches a maximum around the ELA; at lower altitudes higher air temperatures and lower snow amounts reduce the snow cold content (refreezing potential) at the start of the melt season, while at higher altitudes water percolates deep in the firn pack and stored irreducible water continuously refreezes during winter, thereby releasing heat and maintaining a relatively warm firn pack with limited refreezing potential at the start of the melt season. This can be seen in Figure 5b, which shows that the most pronounced refreezing of stored irreducible water occurs in the accumulation zone. This occurs mainly during winter, as the cold wave enters the temperate firn. Some refreezing of irreducible water also occurs at lower altitudes during melt–freeze cycles in the course of the melt season. Refreezing of stored slush water (mean 0.02 m w.e. a⁻¹) is most pronounced in the lower accumulation area (Fig. 5c) and is a minor contribution to total refreezing (0.30 m w.e. a⁻¹). Altogether, refreezing of percolating water (mean 0.13 m w.e. a⁻¹) dominates in the ablation area, whereas stored water refreezing is most prominent in the accumulation zone (mean 0.18 m w.e. a⁻¹).

Fig. 5. Long-term mean spatially distributed patterns of refreezing of percolating water (a), refreezing of irreducible water (b), refreezing of slush water (c) and SI formation (d), averaged over the simulation period 1961–2012.

SI forms when extensive melting, ponding and refreezing produces an ice layer on top of the previous year’s summer surface. SI formation in Svalbard has received considerable attention (e.g. König and others, 2002; Wright and others, 2007; Brandt and others, 2008) as it can comprise a substantial component of the annual mass balance, and because it can be challenging to account for in stake-based mass-balance estimates. We quantify the formation of SI by summing the simulated annual amount of ice formed on top of the previous year’s summer surface, which is exposed at the surface during some time in the melt season. We only account for SI that survives the summer melt season and is associated with remaining winter accumulation, and thereby ignore the temporarily formed SI in the ablation zone. The resulting annual mean pattern of simulated SI in Figure 5d illustrates substantial SI formation in the lower accumulation zones of KNG and HDF. On average, SI is equal to 0.04mw.e.a⁻¹, but local amounts may exceed 0.2 mw.e. a⁻¹. Brandt and others (2008) used GPR to reconstruct SI in the lower accumulation zone on Kongsvegen back to the 1960s and estimated a mean SI contribution of 0.16 ± 0.06 m w.e. a⁻¹, which is slightly higher than the 0.09–0.15 m w.e. a⁻¹ we find in the same region. Simulated small-scale spatial variability in SI can in part be ascribed to variations in the local surface slope, affecting firn water runoff rates. The highest SI values are found on flat slopes, where runoff is limited and ponding water can refreeze during cold periods.

4.1.2. Interannual variability

The area-averaged mass-balance time series (Fig. 6a) is similar for KNG and HDF, and is characterized by high interannual variability. Figure 6b shows the annual time series of spatially averaged annual precipitation and summer (June, July, August (JJA)) air temperature. As in Van Pelt and others (2012), interannual mass-balance variability is strongly anticorrelated with summer temperature (R = −0.73) and positively correlated with annual precipitation (R = −0.61). Accordingly, the more negative mass balances observed since 2000, relative to the preceding decades (1961–99), coincide with a period of decreasing annual precipitation (on average, 90 mm a⁻¹ lower after 2000) and higher summer temperatures (on average, 0.51 K higher after 2000).

Fig. 6. Time series of the area-averaged annual net mass balance for KNG and HDF (a), and the area-averaged annual mean precipitation and summer (JJA) surface temperature (b). The years on the x-axis cover the period from 1 September (preceding year) to 31 August.

Time series of melt and runoff in Figure 7a show that a substantial fraction of melt does not run off as water and is retained instead due to refreezing or (temporary) storage in snow or firn (Fig. 7b). While melting shows strong interannual variability and a clear positive trend, refreezing is more stable and decreases slightly with time.

Fig. 7. Time series of area-averaged annual melt and runoff (a), and annual refreezing and internal accumulation (b). The years on the x-axis cover the period from 1 September (preceding year) to 31 August.

In contrast to refreezing, internal accumulation is found to vary substantially from year to year, with annual mean values ranging from 0.06 to 0.25 m w.e. a⁻¹ (Fig. 7b). In years with little snow accumulation, meltwater is more likely to percolate and refreeze below the previous year’s summer surface, thereby explaining the apparent anticorrelation between annual precipitation (Fig. 6b) and internal accumulation (Fig. 7b). Additionally, cold winter conditions increase the cold wave’s penetration depth and enhance internal accumulation. Thus, internal accumulation maxima occur in years with cold, dry winter conditions. The above results imply that the underestimation of the observed mass balance from stakes in the accumulation zone, as a result of neglecting internal accumulation, is highly variable from year to year.

4.2. Firn evolution

4.2.1. Multi-decadal changes

By comparing output for the period 2000–12 to the preceding period 1961–99, we can detect recent and ongoing changes in surface/firn conditions and discuss the role of feedbacks.

Figure 8 shows the mean distributed fields of the surface albedo, melt and refreezing for the periods 1961–99 and 2000–12. As discussed in Section 4.1, recent enhanced surface mass loss with respect to the preceding decades coincides with higher summer temperatures and reduced precipitation. However, the change in melt rates of 21% between 1961–99 and 2000–12 (Fig. 8c and d) cannot be ascribed wholly to changing external climate conditions. When elevated melt rates persist for multiple years, there follows both a gradual retreat of the firn line in the accumulation zone and a lengthening in the period of bare-ice exposure in the ablation zone. Both effects tend to reduce the mean surface albedo (Fig. 8a and b), thereby increasing absorption of solar radiation and enhancing further surface melt. The above feedback mechanism, commonly referred to as the snow/ice–albedo feedback, effectively enhances melt rates in the ablation area, but is most pronounced in areas which used to be in the accumulation zone and recently became part of the ablation area (Fig. 8a–d). In response to the shrinking of the accumulation area, lower winter snow accumulation, and warmer winter conditions, refreezing is found to decrease by 8% between 1961–99 and 2000–12 (Fig. 8e and f). The drop in refreezing is substantial in the lower ablation area and most pronounced in the former lower accumulation zone. As a result of the increase in melt (21%) and the drop in refreezing (8%), runoff is found to increase by as much as 31% between 1961–99 and 2000–12 (Fig. 8g and h).

Fig. 8. Comparison of the spatially distributed albedo (a, b), surface melt (c, d), refreezing (e, f) and runoff (g, h) averaged over the periods 1961–99 (top) and 2000–12 (bottom).

To further illustrate the multi-decadal firn development, Figure 9 shows the firn density evolution at H5 and H7 on Holtedahlfonna (Fig. 1). At H7, melt–freeze cycles lead to preferential near-surface refreezing and the formation of ‘summer surfaces’, which propagate downwards in time in response to surface accumulation. Substantial surface melt in the early 1970s led to firn thinning at H5 and near-surface densification at H7. In the following decades (till the end of the 1990s) surface melt was limited, and annual precipitation rates were relatively high, implying firn thickening at H5 and reduced near-surface densities at H7. Since 2000, elevated melt rates have led to rapid thinning of the firn, implying complete firn removal at H5, and strong near-surface densification at H7. The suggested formation of a thick near-surface ice layer in the firn at H7 could act as a barrier for vertical water transport and is likely to promote near-surface lateral water flow on top of the ice layer. Effectively, a thick ice layer may limit deep firn water storage.

Fig. 9. The 1961–2012 subsurface density evolution at stakes H5 (a) and H7 (b) along the centre line of Holtedahlfonna (Fig. 1).

Deep (>10 m) firn temperatures are found to be temperate throughout the year for most of the accumulation area, and no clear trends are apparent over the simulation period. Temperate accumulation zones have previously been observed around Kongsfjorden by Björnsson and others (1996), and are common in many polythermal glaciers around Svalbard (Blatter and Hutter, 1990; Pettersson, 2004). Temperate firn conditions in the accumulation zone imply that annual refreezing is controlled by the degree of seasonal winter cooling at the surface. More specifically, winter surface temperatures determine the degree of refreezing of stored firn water, and the cold wave’s penetration depth defines the potential for refreezing of percolating water at the start of the melt season. Deep temperate firn conditions further imply the possibility of deep firn water storage in the form of a perennial firn aquifer. The presence of a dynamic firn aquifer in the upper reaches of HDF has recently been confirmed by Christianson and others (2015).

4.2.2. Initialization sensitivity

As noted, subsurface conditions at the start of the simulation affect the deep firn evolution during the first model decades, which provided the main motivation for performing model initialization (Section 3.5). The state of the firn at the start of the simulation depends ultimately on the chosen climate conditions during this spin-up. To quantify the sensitivity of the model results to the spin-up climate, sensitivity experiments were performed in which the spin-up temperature and precipitation were perturbed. To reduce the numerical cost, sensitivity runs were only done for the stake locations (K1–9, KR1–3 and H1–10; Fig. 1). Results are summarized in Table 2.

Table 2. Sensitivity of the mass balance (m w.e. a⁻¹), refreezing (m w.e. a⁻¹) and firn air content (m) to perturbed spin-up air temperature T and precipitation P

Positively perturbing spin-up temperature by 1 K leads to a smaller firn area at the start of the simulation. The negative impact on the mass balance is somewhat significant after 1 year (−0.05 m w.e. a⁻¹), but is quickly reduced to more modest values after 5 years (−0.02 m w.e. a⁻¹). Refreezing appears to be mostly insensitive to the spin-up climate. On the other hand, the firn air content is found to recover slowly from a perturbed initial state, as significant offsets are apparent even after 20 years into the simulation. Similar sensitivity results, but of opposite sign, are found when positively perturbing spin-up precipitation by 10% (Table 2). Note that the response time of firn conditions to a change in the surface forcing increases with distance to the surface; this explains the quick response time of the mass balance and refreezing processes, since they are determined mainly by near-surface conditions. Because the main focus of this study is on near-surface processes, we conclude that the uncertainty in the spin-up climate has a minor impact on the presented results.

4.2.3. Validation of firn densities

Density profiles from shallow firn cores, drilled in the accumulation zone on Kongsvegen (K8) in spring 1996, 2001 and 2007 and on Holtedahlfonna (H10) in spring 2003 and 2008, were obtained using a PICO corer to depths ranging from 6 to 14 m. In Figure 10, modelled vertical density profiles at observation dates are validated against measured densities. The model appears to do a good job of simulating general density characteristics at both K8 and H10. In addition, interannual variability in mean firn densities seems to be well captured by the model. It is noteworthy that, in a repeated model experiment without refreezing, densities at 5 m depth are, on average, 120 kg m⁻³ lower for the cores at K8 and 128 kg m⁻³ lower for the H10 cores, compared to the original experiment with refreezing. Hence, without refreezing, substantial offsets between simulated and observed densities would become apparent, thereby confirming the substantial impact of refreezing on the subsurface density structure.

Fig. 10. Comparison of modelled and observed density profiles at K8 (a–c) and H10 (d–e).

A known deficiency in the current generation of firn models is their inability to accurately model microscale vertical density variability. To do so would require more detailed accounting for surface–atmosphere interactions (local accumulation variability, wind crust formation), modelling the microstructure development of snow/firn grains and implementing a water transport scheme that accounts for lateral water flow and small-scale processes such as flow impedance and piping (e.g. Pfeffer and Humphrey, 1998; Bell and others, 2008; Wever and others, 2014). Given the regional scale of the experiments, incorporating such processes is beyond the scope of this work, but may be of relevance in future studies focusing in more detail on the local-scale evolution of firn densities.

5. Conclusions

In this study, we applied a coupled surface energy-balance/ firn modelling approach to simulate the long-term (1961–2012) evolution of the surface mass balance and firn evolution of glaciers around Kongsfjorden. Regional climate model output is used to generate a downscaled climate forcing on a 100 m × 100 m grid. Calibration of model parameters is done to increase consistency of the model output and in situ observational data, while extensive spin-up is performed to generate initialized firn conditions at the start of the model simulation.

Results indicate a slightly positive mass balance for both KNG (0.01 m w.e. a⁻¹) and HDF (0.13 m w.e. a⁻¹) for 1961–2012, which for HDF cannot balance mass losses by calving. Refreezing amounts to 0.30 m w.e. a⁻¹, thereby providing a strong buffer for mass loss. Both refreezing of percolating water (0.13 m w.e. a⁻¹) in spring/summer, and refreezing of stored irreducible and slush water (0.18 m w.e. a⁻¹) in autumn/winter, contribute substantially to total refreezing. Internal accumulation of, on average, 0.11 m w.e. a⁻¹ is equivalent to 13% of annual precipitation and leads to a negative bias of mass-balance measurements with stakes in the accumulation zone. Superimposed ice is found to form in the lower accumulation zones of HDF and KNG and contributes, on average, 0.04 m w.e. a⁻¹ to the mass balance, with highest values (up to 0.25 m w.e. a⁻¹) found on flat terrain.

A comparison of the periods 1961–99 and 2000–12 reveals 21% higher melt rates since 2000, which can in part be ascribed to recent warming and lower precipitation rates. Simultaneously, the firn line retreated substantially, which led to surface albedo lowering and further amplification of melt rates by means of the snow–albedo feedback. Refreezing, which is most pronounced in the accumulation zone, is found to be 8% smaller since 2000 in response to shrinking of the firn area and reduced winter cooling. As a result of lower refreezing rates and accelerated melt, runoff increased by as much as 31% since 2000. The above may serve as an example of how feedbacks can act to amplify mass loss in a warming Arctic environment.

By matching simulated and observed quantities, parameter calibration has helped to minimize uncertainty related to inaccuracy of the model physics, model parameter set-up and climate forcing. Sensitivity experiments with perturbed spin-up climate conditions demonstrate that long-term mass balance and refreezing are largely insensitive to uncertainty in the initialization climate. Furthermore, validation of modelled vertical density profiles against shallow core measurements indicates the model’s ability to reasonably replicate spatio-temporal density variations.

The largest uncertainty in the presented results likely comes from a lack of accounting for small-scale variability in precipitation. In mountainous terrain, snow deposition variability can be large since local topography and wind speed and direction exert a strong influence on where snow falls and how it is redistributed by the wind (e.g. Hodgkins and others, 2005; Lehning and others, 2008; Dadic and others, 2010). As discussed in Nuth and others (2012), possible biases in computed mass-balance estimates may arise as centre-line stake measurements may not be representative over the entire glacier surface. Furthermore, ignoring small-scale accumulation variability not only affects the local surface mass balance, but may also induce a bias in the computed area-averaged mass balance, due to the nonlinear dependence of the mass balance on snow accumulation (Van Pelt and others, 2014). Future work will focus on ways to account for the impact of small-scale accumulation variability in mass-balance models. A possibility is to use a snow transport model to generate precipitation fields, accounting for the interaction of local winds and terrain with snowdrift sublimation and deposition.