The climatic mass balance of Svalbard glaciers: a 10-year simulation with a coupled atmosphere–glacier mass balance model

Abstract

In this study we simulate the climatic mass balance of Svalbard glaciers with a coupled atmosphere–glacier model with 3 km grid spacing, from September 2003 to September 2013. We find a mean specific net mass balance of$-$257 mm w.e. yr${}^{-\mathrm{1}}$, corresponding to a mean annual mass loss of about 8.7 Gt, with large interannual variability. Our results are compared with a comprehensive set of mass balance, meteorological, and satellite measurements. Model temperature biases of 0.19 and $-$1.9 ${}^{\circ }$C are found at two glacier automatic weather station sites. Simulated climatic mass balance is mostly within about 100 mm w.e. yr${}^{-\mathrm{1}}$ of stake measurements, and simulated winter accumulation at the Austfonna ice cap shows mean absolute errors of 47 and 67 mm w.e. yr${}^{-\mathrm{1}}$ when compared to radar-derived values for the selected years 2004 and 2006. Comparison of modeled surface height changes from 2003 to 2008, and satellite altimetry reveals good agreement in both mean values and regional differences. The largest deviations from observations are found for winter accumulation at Hansbreen (up to around 1000 mm w.e. yr${}^{-\mathrm{1}}$), a site where sub-grid topography and wind redistribution of snow are important factors. Comparison with simulations using 9 km grid spacing reveal considerable differences on regional and local scales. In addition, 3 km grid spacing allows for a much more detailed comparison with observations than what is possible with 9 km grid spacing. Further decreasing the grid spacing to 1 km appears to be less significant, although in general precipitation amounts increase with resolution. Altogether, the model compares well with observations and offers possibilities for studying glacier climatic mass balance on Svalbard both historically as well as based on climate projections.

Introduction

The Svalbard archipelago has a glacierized area of ca. 34 000 km${}^{\mathrm{2}}$ (Nuth et al., 2013), representing $\sim$ 4 % of the world's land-ice mass outside the Greenland and Antarctic ice sheets. If completely melted, the glaciers on Svalbard could potentially contribute to sea level rise of 17 $\pm$ 2 mm sea level equivalent (Martin-Espanol et al., 2015). The archipelago has already experienced significant warming during the 20th century (Førland et al., 2011) and, with the expected retreat of the sea ice margin, further warming as well as precipitation increases are expected (Day et al., 2012). Projections presented in the latest assessment report of the IPCC (AR5) shows that annual-mean temperatures in this region could rise between 7 and 11 ${}^{\circ }$C by the end of the 21st century under the RPC8.5 scenario, accompanied by a projected precipitation increase between 20 and 50 % (IPCC, 2013). Svalbard glaciers are therefore expected to undergo significant changes during this century (Day et al., 2012; Lang et al., 2015a). However, reliable estimates of future glacier changes require modeling tools that are able to reproduce recent observations. Current model estimates based on global climate data sets (Marzeion et al., 2012, 2015) show significantly more negative mass balance in this region than satellite altimetry and satellite gravimetry over the last decade (Moholdt et al., 2010; Matsuo and Heki, 2013).

Regional model estimates of surface mass balance of Svalbard glaciers have so far mainly focused on individual or a few glaciers and, as noted by Lang et al. (2015b), typically been based on empirical or statistical models. A number of dynamical downscaling simulations focusing on Svalbard glaciers have been performed (Day et al., 2012; Claremar et al., 2012; Lang et al., 2015a, b). However, only two of these studies compare their output with mass balance observations: Day et al. (2012; hereafter DA12) compare precipitation from HadRM3 RCM (25 km grid spacing) with surface mass balance estimates from Pinglot et al. (1999), and Lang et al. (2015b; hereafter LA15b) compare output from the MAR model (10 km) to Pinglot et al. (1999, 2001) as well as a number of altimetry and gravimetry studies of Svalbard glacier mass balance (Wouters et al., 2008; Moholdt et al., 2010; Nuth et al., 2010; Mémin et al., 2011). Both studies show fair agreement with multi-year accumulation records from Pinglot et al. (1999, 2001). LA15b also find mean elevation changes in good agreement with satellite estimates, even though the differences are substantial for some regions. However, DA12 and LA15b do not validate mass balance estimates on timescales shorter than 4 years, nor do they validate on spatial scales that can capture variations on individual glaciers. DA12 also suggest that a grid spacing of 1–5 km may be needed to simulate surface mass balance in the complex terrain that is typical for Svalbard.

In this study, we aim to (i) simulate the climatic mass balance (CMB) of Svalbard glaciers with a higher resolution than previously used with dynamical downscaling, (ii) validate the model with an extensive set of observations in this region, and (iii) investigate the spatial resolution needed to describe the observations. We apply a coupled atmosphere–glacier mass balance model to the entire Svalbard region with a horizontal grid spacing of 3 km, thereby capturing both regional averages for the period 2003–2013 as well as temporal and spatial variations of individual glaciers. The results are validated with (i) observations from weather stations, (ii) mass-balance stakes from four glaciers, (iii) snow accumulation across Austfonna, measured by ground-penetrating radar (GPR), and (iv) satellite altimetry. To examine the importance of model resolution in this region we also compare results from domains with 9 and 3 km grid spacing and, for a selected month, precipitation results from 9, 3, and 1 km grid spacing domains. Through this high-resolution simulation and extensive model evaluation, we aim to provide a detailed estimate of recent CMB of Svalbard glaciers, including its spatial and temporal variations.

Main figure: land areas in the 3 km model domain. Colors indicate glacier grid cells in different subregions and gray indicates non-glacier land grid cells. Stake and AWS locations at the four main validation glaciers are shown as red ${}^{\ast }$ and black ${}^{\ast }$, respectively. NW: northwestern Spitsbergen; NE: northeastern Spitsbergen; SS: southern Spitsbergen; BE: Barentsøya and Edgeøya; VF: Vestfonna; AF: Austfonna. Overlay figure: glacier hypsometry from 3 km model domain compared with 90 m digital elevation model (Nuth et al., 2013).

[Figure omitted. See PDF]

Methods

In the following sections, we describe the two components of the coupled modeling system: the Weather Research and Forecasting (WRF) model (Sect. 2.1) and the CMB model (Sect. 2.2), including optimizations for high-Arctic conditions made in this study. In Sect. 2.3, we describe the different validation data and sites before clarifying comparison methods in Sect. 2.4. Throughout this study we distinguish between the surface mass balance (SMB) and the CMB as recommended by Cogley et al. (2011). The SMB specifies mass changes between the surface and the last summer surface, whereas the CMB also accounts for internal accumulation and ablation (i.e., below the last summer surface). We consider internal ablation as negligible and it is therefore not explicitly treated in our application.

The Weather Research and Forecasting model

The WRF model is a state-of-the-art mesoscale atmospheric model (Skamarock and Klemp, 2008) widely used for research and forecasting applications. In Svalbard, the model has been used to study both atmospheric boundary layer processes (Kilpeläinen et al., 2011, 2012) and atmosphere–land surface interactions over both tundra (Aas et al., 2015) and glaciers (Claremar et al., 2012), with horizontal grid spacing ranging from several tens of kilometers to sub-kilometer scales. In this study, we use the advanced research WRF version 3.6.1 configured with two nested domains of 9 and 3 km horizontal grid spacing. For a single month we also simulate the main regions of interest with additional nested domains at 1 km. The WRF model setup and forcing strategy follows that of Aas et al. (2015). The outer domain (9 km) covers a region of 1080 $\times$ 1080 km. As boundary conditions for this domain we use the ERA-Interim reanalysis (Dee et al., 2011) with a 6 h temporal resolution. We employ the default boundary configuration in WRF, with the outermost grid-point specified, followed by a four-grid-point relaxation zone. From the outer domain the model is one-way nested down to the 3 km domain covering all of Svalbard (Fig. 1). Within both domains the model is allowed to freely evolve (i.e., no nudging or re-initialization), and sea surface temperatures and sea ice fractions are prescribed based on the OSTIA data set (Donlon et al., 2012). The physical parameterization options in WRF follow Aas et al. (2015) with the exceptions of the boundary layer and surface layer parameterizations, the vertical model resolution, and the use of explicit 6th order horizontal advection diffusion, which are selected following Collier et al. (2013). In addition, we use the newer NoahMP land surface scheme to simulate surface conditions and fluxes at non-glaciated grid cells, as it includes improved snow physics and multiple layers in the snowpack over the original Noah scheme (Niu et al., 2011). All three model domains use the same vertical resolution (40 eta layers up to a model top of 25 hPa) and physical parameterizations, with the exception of the cumulus convection scheme that is only employed in the outer (9 km) domain.

The simulation was performed as one transient 10-year simulation, with a 1-day spinup for the atmosphere. The glacier grid points were initialized with output from an earlier simulation with another version of the WRF-CMB model. This did not include results valid for 2003, so a representative year (2007) was used for initialization. Although not accurate, this was considered far more realistic than the default uniform initialization. The 1-month sensitivity simulation with 9, 3, and 1 km grid spacings has been performed as an isolated additional simulation with initialization directly from ERA-Interim to avoid problems with different spinup times for the different domains.

The glacier CMB model

The land surface schemes in WRF have become more advanced in terms of representing snow processes in the recent years but still only simulate a few snow layers (up to three for NoahMP). They therefore have limitations when it comes to simulating the development of deep multiyear snow packs in the accumulation area of Svalbard glaciers, as well as realistically representing different glacier facies (snow, firn, and ice) during the ablation season. We therefore used a modified version of the glacier CMB model of Mölg et al. (2008, 2009) to simulate glacier grid cells. The CMB model uses near-surface temperature, humidity, pressure, and winds, as well as incoming radiation and precipitation as input from WRF, and computes the column specific mass balance from solid precipitation, surface and subsurface melt, refreezing, and liquid water storage in the snowpack, and surface vapor fluxes. The model solves the surface energy balance (SEB) to determine the energy available for surface melt, and resolves the glacier subsurface down to a user defined depth (here 20 m divided into 17 vertical layers). For glacier grid cells, the CMB model updates surface mass and energy fluxes in the WRF model, as well as surface temperature, roughness, and albedo, resulting in a two-way coupled model system (WRF-CMB). Further information about the interactive coupling between the CMB model and WRF is given by Collier et al. (2013, 2015).

We make two adjustments to the CMB model to improve its suitability for Svalbard conditions. First, the albedo scheme (originally based on Oerlemans and Knap, 1998) has been modified to include a separate value for firn albedo (in addition to the standard categories of fresh snow, old snow, and ice). Firn is defined here as snow remaining from the last summer. Secondly, we introduce different aging models for snow (which determines the transition from fresh to old snow) for melting and sub-melting temperatures, by multiplying the snow aging parameter by a factor (warmfact, Table 1) when skin temperatures are at the melting point. The various albedo parameters (Table 1) have been selected based on observations at Austfonna and thereafter adjusted to improve the simulated summer mass balance compared to in situ observations during 1 test year (2005–2006; see also Sect. 6.1.1).

Albedo parameters used in the CMB model.

Parameter	Value
Ice albedo	0.33
Firn albedo (new)	0.50
Old snow albedo	0.63
Fresh snow albedo	0.87
Timescale	15 days
warmfact (new)	5
Depth scale${}^{*}$	3.0 cm

${}^{*}$ Refers to physical snow depth.

Description of observations.

Region	Data	Time period	References
Austfonna, including	AWS	Hourly, since 2004; some data gaps	Schuler et al. (2013)
Etonbreen
	MB stakes	Annually, since 2004	Moholdt et al. (2010); Østby et al. (2013)
	GPR	Annually, since 2004	Taurisano et al. (2007); Dunse et al. (2009)
	Snow pits	Annually, since 2004	Dunse et al. (2009)
Kongsvegen	AWS	Hourly, since 2004; some data gaps	Karner et al. (2013)
	MB stakes	Biannually, since 1986	Nuth et al. (2012)
Holtedahlfonna	MB stakes	Biannually, since 2003	Nuth et al. (2012)
Hansbreen	MB stakes	Biannually, since 1988 (2 years missing).	Grabiec et al. (2012)
		Weekly/monthly since 2005
Svalbard	Satellite altimetry		Moholdt et al. (2010)

Validation data and sites

Table 2 provides an overview of the data sets used for model evaluation. We selected four glaciers, described in more detail below, all of which have mass balance stake measurements at six or more locations, covering all or most years in our study period. The selected sites represent different conditions in terms of glacier geometry, geographical location, local meteorology, altitudinal range, and spatial extent (Fig. 1). Annual measurements of stake heights above the snow surface, snow depth, and density yield specific values of summer and winter-mass balance ($b_{\text{s}}$ and $b_{\text{w}}$ ), which are combined to give the specific net mass balance, $b_{\text{n}}$ (SMB), for each balance year (i.e., between two consecutive end of summers).

GPR measurements of snow accumulation along several transects across the Austfonna ice cap were made each spring in the period 2004–2013. Snow water equivalent values are derived from the radar estimated snow depths multiplied by snow density determined at several snow pits, as described in more detail by Dunse et al. (2009).

Meteorological records at hourly resolution are available from two automatic weather stations (AWSs), one at Austfonna and one on Kongsvegen. Both data sets contain some data gaps. We extract hourly measurements of air temperature ($\sim$ 2 m) and radiation (incoming and outgoing short- and longwave), to calculate daily means, excluding days with incomplete records. Albedo is calculated as the ratio of the outgoing (SW${}_{\mathrm{out}}$) to incoming (SW${}_{\mathrm{in}}$) shortwave radiation, excluding observations outside the range [0.15, 0.95] or with SW${}_{\mathrm{in}}$ < 10 W m${}^{-\mathrm{2}}$. We estimate daily mean albedo using the five hourly observations closest to solar noon to minimize the effect of low solar angle with associated large variations in measured albedo (Schuler et al., 2013).

A geodetic mass-balance estimate from repeat satellite altimetry for the period 2003–2008 (Moholdt et al., 2010) serves as an independent validation of the surface height changes due to climatic mass balance processes. However, it should be kept in mind that the geodetic mass balance reflects both climatic and dynamic mass balance, i.e., it includes mass transfer from higher to lower elevations and losses due to calving at marine termini. We use the regional mean values between 2003 and 2008 according to Fig. 1.

Bias, mean absolute error (MAE), and correlation (corr.) between simulated and observed daily mean temperature, radiation, and albedo at Etonbreen and Kongsvegen.

	Etonbreen			Kongsvegen
	Bias	MAE	Corr.	Bias	MAE	Corr.
T${}_{\mathrm{2}}$ (${}^{\circ }$C)	$-$1.9	2.4	0.96	0.19	1.4	0.98
SW${}_{\mathrm{in}}$ (W m${}^{-\mathrm{2}}$)	$-$13	43	0.81	$-$0.28	21	0.96
SW${}_{\mathrm{out}}$ (W m${}^{-\mathrm{2}}$)	$-$10	17	0.93	$-$6.7	19	0.96
LW${}_{\mathrm{in}}$ (W m${}^{-\mathrm{2}}$)	$-$12	25	0.82	$-$6.9	17	0.89
LW${}_{\mathrm{out}}$ (W m${}^{-\mathrm{2}}$)	$-$3.8	11	0.95	0.32	6.1	0.97
Albedo	$-$0.10	0.12	0.73	$-$0.05	0.08	0.80

Austfonna ice cap, Northeast Svalbard

The Austfonna ice cap in the northeastern part of Svalbard (centered at 79.7${}^{\circ }$ N, 24.0${}^{\circ }$ E) is the largest ice cap of the archipelago. It covers an area of 7800 km${}^{\mathrm{2}}$ and has a simple dome-shaped topography, rising from sea level up to an elevation of $\sim$ 800 m a.s.l. (Moholdt and Kääb, 2012). The recent CMB of Austfonna was nearly in balance (Moholdt et al., 2010), yet the ice cap was losing mass due to calving and retreat of the marine margin (Dowdeswell et al., 2008). Snow accumulation is spatially and temporally heterogeneous and asymmetrical across the ice cap, with amounts in the southeast being double the amounts in the northwest and large interannual variability along all profiles (Pinglot et al., 1999; Taurisano et al., 2007; Dunse et al., 2009).

Since spring 2004, field measurements have been performed annually by the University of Oslo and the Norwegian Polar Institute. Available data include records from about 20 mass balance stakes, annually repeated GPR and kinematic GNSS (Global Navigation Satellite System) profiling across the ice cap, and snow pit investigations of snow depth and density (Taurisano et al., 2007; Dunse et al., 2009). In the present study we compare GPR-derived winter accumulation with the corresponding WRF-CMB results, averaging all available in situ measurements within a particular WRF-CMB grid cell (Sect. 3.3).

Etonbreen is located at the western part of Austfonna, with six stakes, and an AWS operated since 2004. The AWS is located at 22${}^{\circ }$25${}^{\prime }$12${}^{\prime \prime }$ E, 79${}^{\circ }$43${}^{\prime }$48${}^{\prime \prime }$ N and 370 m a.s.l., just below the mean equilibrium line altitude (ELA) of $\sim$ 400 m (Schuler et al., 2013).

Simulated (blue) and observed (red) albedo at Etonbreen and Kongsvegen for the summers 2008 and 2013, with simulated solid (blue) and liquid (yellow) precipitation indicated as bars.

[Figure omitted. See PDF]

Kongsvegen and Holtedahlfonna, northwestern Spitsbergen

Kongsvegen (78.8${}^{\circ }$ N, 13.0${}^{\circ }$ E) and Holtedahlfonna (79.0${}^{\circ }$ N, 13.5${}^{\circ }$ E) are both located near Ny-Ålesund, in northwestern Spitsbergen.

Kongsvegen is a $\sim$ 100 km${}^{\mathrm{2}}$, $\sim$ 27 km long valley glacier extending from an ice divide at $\sim$ 800 m a.s.l. down to sea level. Outflow at its marine terminus is restricted by its fast-flowing neighbor Kronebreen, with which Kongsvegen shares a small fraction of the calving front. The Norwegian Polar Institute has measured winter and summer mass balance at nine stakes since 1986 (Hagen et al., 2003; Nuth et al., 2012; Karner et al., 2013). Kongsvegen is a surge-type glacier, currently in its quiescent phase, since the last surge around 1948. Observed elevation changes are dominated by the SMB (Melvold and Hagen, 1998).

North of Kongsvegen is Holtedahlfonna, the upper catchment of the Holtedahlfonna–Kronebreen glacier system, whose total area is $\sim$ 390 km${}^{\mathrm{2}}$, and which extends up to an elevation of $\sim$ 1400 m a.s.l. SMB has been studied on Holtedahlfonna since spring 2003, using 10 mass-balance stakes.

Stakes at both glaciers are measured in nearly all years in spring and at the end of summer, typically in early September. During the last 2 decades, the SMB of both glaciers has changed from close to 0 to increasingly negative values (Kohler et al., 2007; Nuth et al., 2012). Since spring 2000, an AWS has been operated on Kongsvegen, at 78.76${}^{\circ }$ N, 13.16${}^{\circ }$ E at 537 m a.s.l., close to the ELA (Karner et al., 2013).

Hansbreen, southern Spitsbergen

Hansbreen (77.1${}^{\circ }$ N, 15.6${}^{\circ }$ E) is located in the southern part of Spitsbergen and covers an area of $\sim$ 56 km${}^{\mathrm{2}}$. It is a 15 km long valley glacier, extending from its $\sim$ 1.5 km wide active calving front in Isbjornhamna, Hornsund, up to an ice divide at $\sim$ 490 m a.s.l. The Hansbreen system consists of the main trunk glacier and four tributary glaciers on the west side (Grabiec et al., 2011). SMB on Hansbreen has been studied since 1989 when 11 mass balance stakes were deployed along centerline of the glacier. Since 2005 these stakes are measured on a weekly basis in ablation zone and on a monthly basis in accumulation area (Grabiec et al., 2012).

Snow accumulation on Hansbreen is highly variable. Earlier studies show that there is a strong asymmetry between eastern and western side of the glacier. Minimal snow accumulation is observed in southern part of the tongue and along eastern side of the glacier. On the eastern side, the glacier is bordered by the massif of Sofiekammen, which forms an orographic barrier for advecting air masses. Therefore, a foehn effect is widely observed during strong easterly winds, causing deflation of snow and redeposition towards the western side of the glacier (Grabiec et al., 2006).

Comparison methods

We compare our model results with AWS (Sect. 3.1) and stake (Sect. 3.2) data using the WRF-CMB grid point nearest to each data point. All stakes on a glacier are compared with the corresponding grid cells to yield information both about mean values as well as gradients along the glacier. The ELA (also Sect. 3.2) is calculated by linear interpolation of the two stakes or grid cells with the least positive and negative mass balance. When all stakes or grid cells have the same sign, we use the maximum and minimum value to extrapolate to the ELA. Note that the CMB simulated by WRF-CMB also includes internal refreezing below the last summer surface, which is not captured by the stake SMB.

Model validation

In the following sections, we assess the model performance by comparing the simulated SEB and temperature with the AWS data (Sect. 3.1), the CMB at the four glaciers with stake measurements (Sect. 3.2), the winter accumulation at Austfonna with GPR measurements (Sect. 3.3), and, finally, the simulated regional height changes over the first 5 years of the study period with the altimetry data (Sect. 3.4). All model results in this section are from the 3 km domain.

Weather stations

The AWSs provide key information about the conditions at the glacier surface throughout the year and therefore permit a detailed evaluation of important aspects of the model. Table 3 compares simulated and observed air temperature, radiation fluxes, and albedo at Kongsvegen and Etonbreen. Note that the model grid points are 44 and 29 m lower than the station heights at these two locations, respectively. At Kongsvegen the simulated temperature compares very well with observations, with a bias of 0.19 ${}^{\circ }$C and correlation of daily mean values of 0.98. At Etonbreen the simulated temperatures are somewhat too low ($-$1.9 ${}^{\circ }$C) but have a similar correlation (0.96). The radiation components also show better agreement at Kongsvegen than at Etonbreen. At Kongsvegen, biases ranging from 0.3 W m${}^{-\mathrm{2}}$ (outgoing longwave, LW${}_{\mathrm{out}}$) to $-$6.9 W m${}^{-\mathrm{2}}$ (incoming longwave, LW${}_{\mathrm{in}}$), whereas the radiation biases at Etonbreen vary from $-$3.8 (LW${}_{\mathrm{out}}$) to $-$13 W m${}^{-\mathrm{2}}$ (SW${}_{\mathrm{in}}$). There is also a noticeable albedo bias of $-$0.10 at Etonbreen, compared to only $-$0.05 at Kongsvegen.

Simulated (blue) and observed (red) annual, winter, and summer mass balance at the four main validation glaciers. The 10-year mean is indicated by solid lines and the years 2004 and 2008 (negative and positive anomaly, respectively) are shown as dashed lines. Stars indicate elevation of individual stakes (red) or grid cells (blue).

[Figure omitted. See PDF]

Estimated ELA from WRF-CMB (blue) and stakes (red) at the four main validation glaciers.

[Figure omitted. See PDF]

(a) Simulated and GPR-derived winter accumulation at Austfonna 1 May 2004. (b) Scatterplot of data from GPR locations in (a). Red dots show points located northwest of summit, and black dots show points located southeast of summit. (c) Same as (a) for 1 May 2006. (d) Same as (b) for 1 May 2006.

[Figure omitted. See PDF]

A more detailed comparison of simulated and observed albedo from the summers of 2008 and 2013 is shown in Fig. 2, along with simulated solid and liquid precipitation. Overall, the model simulates well both the magnitude and temporal changes in albedo. The close connection between simulated solid precipitation events and observed albedo increase indicates that the model adequately captures both the timing and phase of the summer precipitation. The simulated albedo response to solid precipitation is also similar to measurements, except at Etonbreen in 2013. Here the model does not simulate a large enough increase in albedo after snow events, which might indicate that the model is not sensitive enough to small amounts of snow on ice (i.e., too large depth scale in Table 1) or underestimates precipitation or its frozen fraction on these occasions. Additionally, the model underestimates snow albedo during much of the summer in 2008, and at Etonbreen also simulates some periods with blue ice. While the observations show a very slow decrease in albedo over this summer, the modeled albedo quickly drops to values below 0.7. This could be related to the almost complete lack of rain events during the summer of 2008. In 2013, however, the model simulates numerous and large rain events during the summer, which seems to coincide with observed drops in albedo. This suggests that snow wetness needs to be accounted for in the albedo parameterization to realistically simulate snow albedo on Svalbard and that the snow aging parameters used here (Table 1) are more appropriate for wet rather than dry summer conditions.

Altogether, the model seems to reasonably reproduce the local conditions at the AWS stations. Individual radiation biases are found, which are likely to impact the quality of the CMB estimations negatively. However, this is a known challenge in the Arctic, where atmospheric models often show significant biases in key cloud properties (Morrison et al., 2009).

Stakes

Figure 3 shows mean annual ($b_{\text{n}}$), winter ($b_{\text{w}}$), and summer ($b_{\text{s}}$) CMB at each stake over the study period. Overall, the model reproduces the mean observed mass balances well, including differences in mean values between the four glaciers and the gradients at each glacier. The main exception is the winter balance at Hansbreen, which is considerably underestimated by the model (see Sect. 4). Hansbreen also shows large variations along the glacier, both in $b_{\text{s}}$ and $b_{\text{w}}$, which the model is not able to reproduce.

Looking at individual years, the model correctly identifies the years 2004 and 2008 as having anomalously low and high mass balance, respectively (Fig. 3, stippled lines). The agreement, however, is not as good for these 2 years as for the mean values.

To look further at temporal variations, we compare modeled ELA to that derived from stake measurements (Fig. 4). In general there is good agreement in terms of both ELA and interannual variability, with the exception again being Hansbreen. Here the model strongly overestimates ELA, in accordance with the underestimation of $b_{\text{w}}$ (Fig. 3). However, at this glacier, the observed mass balance–elevation relationship shows considerable nonlinearity, with the observations indicating 0 $b_{\text{n}}$ at several elevations during some years, rendering ambiguous ELA estimates. For the other three glaciers the model simulates ELA well, including the large difference between Kongsvegen and Holtedahlfonna, which are located close to each other (Fig. 1).

Mean surface elevation change rates (m yr${}^{-\mathrm{1}}$) for 2003–2008 from WRF-CMB and Moholdt et al. (2010). The Svalbard mean value from Moholdt et al. (2010) used here is the mean of the regions included here (Fig. 1) weighted by area.

	WRF-CMB	Moholdt et al. (2010)
Northwestern Spitsbergen	$-$0.44	$-$0.54
Northeast Spitsbergen	0.02	0.06
Southern Spitsbergen	$-$0.33	$-$0.15
Barenstøya/Edgeøya	$-$0.35	$-$0.17
Vestfonna	$-$0.34	$-$0.16
Austfonna	0.05	0.11
Svalbard	$-$0.18	$-$0.11

GPR-derived snow profiles

To evaluate the winter-mass balance simulated by WRF-CMB, we compare simulated precipitation across Austfonna between early September and early May with observations of snow accumulation by GPR. The results from May 2004 and May 2006 (Fig. 5) demonstrate the large interannual variability in the total amount of snow, with 2004 and 2006 representing a low (Fig. 5a) and high (Fig. 5c) accumulation year, respectively. WRF-CMB captures the spatial pattern and total amount of snow very well, with average biases of about 2 and 6 % and mean absolute errors (MAEs) of 47 and 67 mm w.e. yr${}^{-\mathrm{1}}$ in 2004 and 2006, respectively. Local differences, for example an overestimation of snow within the ablation area towards the lower end of the western profile in both years (Fig. 5b, d), may partly be explained by wind erosion, which is not represented in the model. The model also tends to slightly underestimate snow accumulation in the summit area, which again could be related to wind redistribution of snow or, by a negative elevation bias of the WRF-CMB digital elevation model which is on average $\sim$ 20–50 m lower than the average position of the GPR measurements.

Mean annual CMB (mm w.e. yr${}^{-\mathrm{1}}$) from 2003 to 2013 simulated with 9 km (left) and 3 km (right) grid spacings.

[Figure omitted. See PDF]

Model topography at 9, 3, and 1 km grid spacings around Etonbreen (top), Kongsvegen and Holtedahlfonna (middle), and Hansbreen (bottom). Red dots indicate stake locations.

[Figure omitted. See PDF]

Precipitation during October 2007 at stake locations with 9 km (blue), 3 km (green), and 1 km (red) grid spacing.

[Figure omitted. See PDF]

Satellite altimetry

We perform a regional evaluation of the simulations by comparing surface height changes with those measured by satellite altimetry in the period 2003–2008 (Moholdt et al., 2010). Note that measured height changes are also a result of glacier dynamics (although the retreat or advance of the calving front was not accounted for). The model, however, does not include any glacier dynamics and the numbers are therefore not directly comparable. Still, both model and satellite data show Austfonna and northeastern Spitsbergen to be the only regions with positive surface height change during these years, and northwestern Spitsbergen as the region with the largest surface lowering (Table 4). The other three regions all show moderate lowering in both estimates. The model therefore seems to capture regional differences very well during this period. The mass loss from calving flux has been estimated to be 6.75 km${}^{\mathrm{3}}$ yr${}^{-\mathrm{1}}$ (w.e.) for the years 2000–2006 (Blaszczyk et al., 2009), which corresponds to an additional lowering of about 0.2 m yr${}^{-\mathrm{1}}$. This suggests that the model in general simulates too much surface lowering in this period. However, the time periods are different and the estimate of Moholdt et al. (2010) does not include the effect of retreat or advance of the calving front. One must therefore still be cautious when comparing these numbers.

Together, these results show that the model captures the spatial and temporal variability in CMB across Svalbard well, with mean values in fair agreement with observations. The largest model errors are found at Hansbreen, which has both the largest cross-glacier variability in accumulation and as a relatively steep accumulation gradient, such that the stake measurements may not correlate as well to the model as at other sites with a greater degree of spatial homogeneity. This raises the question about the sensitivity of model resolution in relation to capturing the main topographic features.

Sensitivity to model resolution

Svalbard topography is relatively rugged, with fjords, tundra, glaciers, and mountains all found in close proximity to each other. Dynamical downscaling requires a discrete representation of the topography and has therefore limited spatial resolution, which potentially can be insufficient to resolve a number of small-scale processes. In this section, we investigate the sensitivity of simulated CMB to model resolution by comparing results from the 9 and 3 km domains for the entire model period. We then evaluate precipitation amounts and distribution in a separate simulation with 9, 3, and 1 km grid spacings (hereafter R9, R3, and R1) of the month of October 2007, which was among the wettest months in the 10-year period, with a number of smaller and larger precipitation events.

Figure 6 shows the mean annual CMB in R9 and R3 over the entire period. When averaged across all of Svalbard, the difference is relatively small ($\sim$ 30 mm w.e. yr${}^{-\mathrm{1}}$). On a regional scale, however, differences are more substantial, with the southeastern islands (BE region in Fig. 1) having almost 100 mm w.e. yr${}^{-\mathrm{1}}$ more negative mass balance in R9 than in R3. The modest difference found when averaging across the entire archipelago is therefore at least partly due to compensating regional differences. On the local scale, differences are even more pronounced. The lack of small glaciers in R9 gives large ice-free areas compared with R3 (e.g., in much of central Spitsbergen). The CMB gradient is also often larger in R3 than R9, as higher maximum and lower minimum values are resolved, which is likely key for long-term model simulations, where geometry adjustments of glaciers are considered.

Mean annual CMB fluxes (mm w.e. yr${}^{-\mathrm{1}}$) and standard deviations (SD) for Svalbard and each subregion in Fig. 1. SPR: solid precipitation; REF: refreeze; SUI: superimposed ice; DEP: deposition; MLT: surface melt; SUM: subsurface melt; SLB: sublimation; DLW: change in snow liquid water.

	Svalbard		NW		NE		SS		BE		VF		AF
	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Mean	SD
SPR	639	107	588	107	647	118	774	138	703	113	491	91	586	111
REF	141	21	126	23	151	23	143	21	125	18	112	23	152	27
SUI	60	9.2	54	10	65	10	61	9.0	54	7.7	48	10	65	12
DEP	51	7.9	50	7.1	45	6.5	64	11	64	13	47	8.7	44	6.8
MLT	$-$1000	262	$-$1090	243	$-$837	250	$-$1350	315	$-$1260	398	$-$876	271	$-$776	256
SUM	$-$124	27	$-$147	38	$-$111	26	$-$148	33	$-$132	36	$-$114	32	$-$98	26
SBL	$-$26	2.6	$-$28	3.3	$-$27	2.3	$-$23	2.6	$-$22	3.1	$-$26	3.2	$-$28	3.2
DLW	4.0	12	3.1	10.3	5.8	16.8	2.2	10.7	1	6.4	2.7	10.2	5.5	19
CMB	$-$257	358	$-$439	376	$-$62	347	$-$481	464	$-$465	484	$-$314	337	$-$49	323

(a) annual mean mass balance (m w.e. yr${}^{-\mathrm{1}}$) for Svalbard (black, left $y$ axis) and regional deviations (colors, right $y$ axis). (b) Same as (a) but for winter-mass balance (September to April). (c) Same as for (a) but for summer mass balance (May to August).

[Figure omitted. See PDF]

The level of detail resolved also plays an important role in model evaluation, especially when comparing to in situ point measurements. Figure 7 shows the terrain height in the three regions with stake measurements for R9, R3, and R1. Here we see new topographic features and a more detailed coastline emerging with each increase in resolution. In R9 the topographic gradients along the glaciers with stakes are mostly too low, and many stakes maps on to the same grid cell. Conversely, in R3 most stakes maps on to individual grid cells and the altitudinal range is in better agreement with the stakes (see also Fig. 3). Still, there are distinct topographic features around these four glaciers that only emerge in R1. For example, the model represents Hansbreen, Kongsvegen, and the upper part of Holtedahlfonna as valley glaciers only at this resolution. The simulated precipitation from October 2007 (Fig. 8) also shows distinct differences between the three resolutions. At Etonbreen, the shape of the precipitation curve along the glacier is very different in R9 compared with R3 and R1. The other three glaciers show more linear increases in precipitation with altitude during this month at all three resolutions. Still the precipitation amount varies and mainly increases with higher resolution, with local differences of up to about 50 mm (25 %) at Hansbreen ($\sim$ 300 m a.s.l.).

We therefore conclude that increasing resolution from 9 to 3 km grid spacings is important for simulating precipitation and CMB on local and regional scales, as well as for reliable comparison with in situ point measurements. Increasing the resolution further to 1 km grid spacing seems to have a smaller effect on the four glaciers investigated, although the precipitation amount is further increased. Model resolution might therefore explain some of the negative model bias in $b_{\text{w}}$ (cf. Fig. 3). However, resolution alone does not explain the large $b_{\text{w}}$ bias at Hansbreen. Instead it is likely that processes like wind drift and snow redistribution are important for the accumulation pattern on this glacier. These processes likely also affect the accumulation at the other glaciers; however, their influence is largely unconstrained.

Climatic mass balance

We now turn to simulated annual and seasonal CMB results for Svalbard and the different subregions (Fig. 1). Consistent with previous studies (e.g., Hagen et al., 2003; Lang et al., 2015b), we find large interannual variation in CMB (Fig. 9 and Table 5), which precludes trend analysis on the timescales considered in this study. Likewise, our mean mass balance value of $-$257 mm w.e. yr${}^{-\mathrm{1}}$ depends on the model period considered; 5-year mean values vary from $-$313 mm w.e. yr${}^{-\mathrm{1}}$ (2009–2013) to $-$27 mm w.e. yr${}^{-\mathrm{1}}$ (2006–2010). Multiplied with a glacier area of 34 000 km${}^{\mathrm{2}}$ (Nuth et al., 2013), these values correspond to a mean annual mass loss of 8.7, 11, or 0.92 Gt, respectively.

Comparing the two seasons reveals that summer mass balance varies about twice as much as winter-mass balance (Fig. 9), indicating that ablation processes vary more from year to year than winter accumulation. This is confirmed from the annual mass fluxes shown in Fig. 10b, with solid precipitation showing much less variability than surface melting. Melting is in turn largely a result of the radiation imbalance during the summer months (JJA; Fig. 10a), which on average accounts for about 80 % of the net energy to the surface during these months. However, years with anomalously large melting (especially 2004 and 2013) are characterized by much larger than normal sensible and latent heat fluxes at the surface. As these fluxes indicate warmer and moister air in the atmosphere relative to the glacier surface, the melting anomalies likely result from advection of warm air from outside the region, as was also found by Lang et al. (2015b) for the year 2013.

For the individual regions we find large differences in simulated CMB that persists throughout the whole period (Fig. 9; right axis). Most noticeably the northeastern regions (AF, VF, and NE) show less negative summer balance than the Svalbard average and mostly lower than average winter accumulation. The regions with high winter accumulation in the south and east correspond to the regions with the lowest ELA reported by Hagen et al. (2003; see also mean annual precipitation and ELA estimates in the Supplement).

(a) Mean summer (JJA) surface energy balance fluxes. $Q_{\text{R}}$: net radiation; $Q_{\text{M}}$: melt energy; $Q_{\mathrm{SL}}$: sensible and latent heat flux; $Q_{\mathrm{PRC}}$: heat from precipitation; $Q_{\text{C}}$: ice heat flux; $Q_{\mathrm{PS}}$: penetrating solar radiation. (b) Annual mass fluxes averaged over Svalbard. SPR: solid precipitation; REF: refreeze; SUI: superimposed ice; DEP: deposition; MLT: surface melt; SUM: subsurface melt; SLB: sublimation; DLW: change in snow liquid water. The resulting CMB is indicated by white dots.

[Figure omitted. See PDF]

Discussion

Comparison with earlier studies

We will focus our comparison here with the two regional model studies that report comparison with SMB or CMB measurements, namely DA12 and LA15b (Sect. 1). Both of these studies compare their results with ice core measurements in accumulation areas (Pinglot et al., 1999, 2001), although DA12 only include accumulation and not melting in their simulation. DA12 find biases ranging from $-$240 to 130 mm w.e. yr${}^{-\mathrm{1}}$ with an MAE close to 100 mm w.e. yr${}^{-\mathrm{1}}$, whereas LA15b report biases between $-$310 and 140 mm w.e. yr${}^{-\mathrm{1}}$, with a MAE of 100 mm w.e. yr${}^{-\mathrm{1}}$. Our simulation period does not cover these ice core measurements so that a direct comparison is not possible. Still, the MAEs in winter accumulation at Austfonna (Fig. 5) are 47 and 67 mm w.e. yr${}^{-\mathrm{1}}$ for 2004 and 2006, respectively. At the stake locations, the simulated CMB (Fig. 3) is mostly within about 100 mm w.e. yr${}^{-\mathrm{1}}$ of the observations, except at Hansbreen where it ranges from close to 0 to about 1000 mm w.e. yr${}^{-\mathrm{1}}$. Our CMB results are therefore an improvement compared to DA12 and LA15b when Hansbreen is ignored (where processes not represented in our model are believed to be important).

The near-surface temperature biases are in both DA12 and LA15b larger than ours. DA12 report annual biases “less than 2 ${}^{\circ }$C” for two stations and “less than 6 ${}^{\circ }$C” for the third (Kongsvegen) and LA15 show annual biases between $-$1.3 and $-$4.0 ${}^{\circ }$C (all negative). We only use data from the two AWSs on the glaciers here, with biases of 0.19 and $-$1.9 ${}^{\circ }$C at Kongsvegen and Etonbreen, respectively. These biases span the 2-year mean biases found by Claremar et al. (2012) and are also similar to those found by Aas et al. (2015) for the year 2008–2009 (including also tundra sites), both using the WRF model.

Some of the improvement found in the present study can probably be attributed to the increased model resolution. Both smaller elevation differences between stations and the model grid cell and better resolved surrounding topography likely improved the results. We also note here that our results do not cover the same time periods as DA12 and LA15b, so that the quality of the boundary conditions might differ. However, LA15b report relatively small biases for ERA-Interim (between $-$1.95 and 2.24 ${}^{\circ }$C) despite similar or larger elevation differences and covering the same period as LA15b. Also, the large increase in resolution from DA12 (25 km) to LA15b (10 km) is not accompanied with a clear improvement in the mass balance simulation. It therefore seems clear that the WRF-CMB model with the setup used here offers a real improvement over DA12 and LA15b.

Model uncertainties

In the following section we discuss the main uncertainties in the model results, starting with the atmospheric forcing, before discussing the representation of albedo, turbulent fluxes, and subsurface processes in the CMB model. Although we recognize that the observations have uncertainties and limitations, it is beyond the scope of this study to go into those here.

Atmospheric forcing

The quality of any CMB simulation depends largely on the atmospheric forcing used. In this study we have used a coupled atmosphere–glacier model in which the atmospheric component (WRF) has been well tested for this region. It has been shown to produce temperatures that are in good agreement with observations and relatively small biases in energy fluxes on annual timescales, although they can differ significantly on seasonal and shorter timescales (Aas et al., 2015). Deviations from observations can come both from insufficient representation of processes within the WRF model and from errors in the boundary conditions (here ERA-Interim and OSTIA). However, comparison with weather stations (Sect. 3.1) and the results from Aas et al. (2015) show that there is good agreement with observations for the atmospheric part compared to similar studies (Claremar et al., 2012; DA12; LA15b).

Albedo

The albedo parameterization in the CMB model uses a set of parameters that likely varies considerably for different locations but is treated as spatially invariant. Initially, these parameters were set based on data from the AWS at Austfonna, where we had the longest series of radiation data. This, however, gave too little summer ablation in general, and these parameters were therefore instead tuned to give summer mass balance values in line with observations based on a set of sensitivity simulations of the year 2005–2006 with the 9 km grid spacing domain. The resulting values (Table 1) are similar to those used by van Pelt et al. (2012). With the complex model system used here we cannot perform a set of simulations at the full resolution and time period, but we can acknowledge that these values are uncertain and should ideally vary across the region. Greuell et al. (2007) found MODIS ice albedo values at different glaciers across Svalbard between 0.44 (Etonbreen) and 0.13 (Lisbethbreen), which also explains why our initial values from Etonbreen gave too little summer ablation in general. A large step forward for CMB modeling of this region would therefore be to include spatially varying albedo parameters based on satellite measurements. Including snow wetness could also offer improvements for simulation of changes in snow albedo with time (see Sect. 3.1).

Atmospheric stability

We have seen that the sensible and latent heat fluxes are an important part of the SEB during the ablation season (Sect. 5). However, the simulation of these fluxes in stable conditions is a major challenge subject to ongoing research (e.g., Holtslag et al., 2013). In our setup of the CMB model, we use a stability correction based on the bulk Richardson number (Braithwaite, 1995). Conway and Cullen (2013) found that this correction gives too low fluxes during stable conditions and low wind speeds at a New Zealand Southern Alps glacier, with simulated turbulent fluxes close to 0 when observations showed values between 50 and 100 W m${}^{-\mathrm{2}}$. Based on their results, and to avoid runaway cooling of the glacier surface during stable conditions in the winter, we limited the reduction in the turbulent fluxes in stable conditions to 30 %, consistent with previous studies (Martin and Lejeune, 1998; Giesen et al., 2009; Collier et al., 2015). Due to these issues, and since they are not measured at the study glaciers, these fluxes also contribute to the overall model uncertainty.

Subsurface processes

As can be seen from Fig. 10b, the subsurface components (refreezing, superimposed ice, subsurface melting, and change in liquid water content) contribute considerably to the total simulated CMB. An important simplification in the CMB model is the use of a bulk snow density, as it is not intended for detailed snowpack studies. For the highest areas on Austfonna and northeastern Spitsbergen, where the simulated snow and firn accumulation reaches values of 14 m during the simulation, this simplification is likely inappropriate. In addition, the CMB model currently calculates superimposed ice as a specified fraction of the internal refreezing of liquid water. Modeling this process is challenging, even with a vertically resolved treatment of snow density. Thus, while the inclusion of the subsurface processes likely offers an important step forward compared to only simulating the surface mass balance, the accuracy of these subsurface processes in the model needs to be improved.

Conclusions

In this study, we simulated the CMB of the Svalbard archipelago with a coupled atmosphere–glacier model, for the period 2003 to 2013 with 9 km and 3 km grid spacings, as well as with 1 km grid spacing for a shorter period. The results have been compared with extensive observational data from weather stations, mass balance stakes, ground penetrating radar and satellite altimetry. Our main findings are as follows.

• The WRF-CMB model with 3 km grid spacing and the configuration used here reproduces observed CMB on Svalbard very well, from local to regional scales.

• We confirm the need for high spatial resolution (1–5 km grid spacing) to realistically simulate CMB at glacier scale, as suggested by Day et al. (2012). The results from 3 and 1 km grid spacings show distinct features at local scales that are not present at 9 km, and mean CMB at regional scale differed by up to $\sim$ 100 mm w.e. yr${}^{-\mathrm{1}}$ between the 3 and 9 km simulations.

• Large variations in CMB on small spatial scales reduce the representativeness of individual point measurements when compared with grid cells larger than 1–5 km.

• We find large year-to-year variability in average CMB on Svalbard during our simulation period, which can be mainly attributed to variations in melting. The largest component in the summer surface energy balance driving this melting is the radiation imbalance, even though temperature-dependent latent and sensible heat fluxes also contribute to much of the year-to-year variability, especially during the years with anomalously large mass loss. More research is, however, needed to better understand the drivers of this variability.

For the first time, this study presents results from dynamical downscaling of Svalbard CMB at the resolution suggested by Day et al. (2012). In addition, we have utilized a large number of observations on different spatial scales to get a robust evaluation of model performance, thereby representing a considerable step forward in the pursuit of reliable simulations of the CMB on Svalbard. Further improvements to several aspects of the WRF-CMB would be desirable for this region, including using spatially variable albedo parameters and improved representation of subsurface processes. Still, the WRF-CMB model – which has not previously been validated for Arctic conditions – has here been shown to be an appropriate tool for studying Svalbard CMB.

The Supplement related to this article is available online at doi:10.5194/tc-10-1089-2016-supplement.

Acknowledgements

We thank two anonymous reviewers for valuable comments and suggestions to the original manuscript and V. Radić for acting as Editor. This study was mainly carried out as a part of the CryoMet project funded by the Norwegian Research Council (NFR 214465). The preparation of this paper has further been supported by the Nordforsk-funded projects Green Growth Based on Marine Resources: Ecological and Sociological Economic Constraints (GreenMAR) and Nordic Center of Excellence sSTICC (eScience Tools for Investigating Climate Change at high northern latitudes; grant 57001), as well as within the statutory activities No 3841/E-41/S/2015 of the Ministry of Science and Higher Education of Poland.Edited by: V. Radic