The mesoscale atmospheric model WRF is used over three Svalbard glaciers. The simulations are done with a setup of the model corresponding to the state-of-the-art model for polar conditions, Polar WRF, and it was validated using surface observations. The ERA-Interim reanalysis was used for boundary forcing and the model was used with three nested smaller domains, 24 and 8 km, and 2.7 km resolution. The model was used for a two-year period as well as for a more detailed study using 3 summer and winter months. In addition sensitivity tests using finer horizontal and vertical resolution in the boundary layer and using different physics schemes were performed. Temperature and incoming short- and long-wave radiation were skillfully simulated, with lower agreement between measured and modelled wind speed. Increased vertical resolution improved the frequency distributions of the wind speed and the temperature. The choice of different physics schemes only slightly changed the model results. The polar-optimized microphysics scheme outperformed a slightly simpler microphysics scheme, but the two alternative and more sophisticated PBL schemes improved the model score. A PBL scheme developed for very stable stratifications (QNSE) proved to be better in the winter.

1. Introduction

The Svalbard archipelago experienced a significant thinning of its glacial mass over the last century, which contributed about 0.026 mm yr−1 sea-level rise over the last 50 years [14]. In consequence, the stratification and circulation of fjord systems respond to enhanced fresh water supply by melt water and icebergs from calving glaciers [5]. Significant changes are also observed regarding Svalbard sea ice, permafrost, or land and marine ecosystems [610]. The thinning and retreat of glaciers on Svalbard is probably both an effect of the warming that set off after the little ice age and the warming of the Arctic over the last decade that likely is an effect of anthropogenic greenhouse gas forcing [11]. The Svalbard archipelago has also experienced a significantly increased precipitation since the 1960s [12]. Several projects have recently been initiated to better understand these detected environmental changes and the response of the Arctic cryosphere (e.g., SVALI; http://www.ncoe-svali.org or SvalGlac; http://svalglac.eu). SVALI (Stability and Variations of Arctic Land Ice) is a Nordic centre of excellence, initiated by the Nordic council of ministers, with the focus to study the interaction between climate change and the cryosphere, with focus on glaciers in the Arctic/North Atlantic area. SVALI’s task is to improve the knowledge of the boundary effects and the processes and their parameterisation between the atmosphere and the cryosphere to deliver better prognostics of the future of the glacial systems in the Arctic.

Precipitation, temperature, humidity, wind, radiation, and turbulent heat fluxes are key meteorological factors determining the mass and energy balance of snow and ice in general. The Arctic regions, however, are prone to a notorious lack of corresponding measurement data. Therefore, GCM and global reanalyses data were early considered as important sources for analysis of weather phenomena and climate in these regions. Meanwhile such is increasingly based on statistically or dynamically downscaled output from regional-scale atmospheric models (RCMs). No doubt, however, observational data still play an important role for process studies, data assimilation, or validation of model results.

The skill of prominent global or regional climate models to simulate the Arctic atmosphere was investigated by numerous validation and model intercomparison studies [1320]. Major efforts were associated to carefully designed and coordinated observational and model intercomparison projects (SHEBA, Surface Heat Budget of the Arctic Ocean, Perovich et al. [21]; ARCMIP, Curry and Lynch [15]). These studies proved the progressive skill of RCMs to simulate the spatiotemporal structure and evolution of the near-surface air properties while revealing some still remaining conflicts on the other hand. Roughly summarizing the intermodel comparisons by, for example, Tjernström et al. 2004 or Rinke et al. 2005 [16, 18], the seasonal weather patterns are skillfully reproduced by the considered models. The simulated near-surface variables (pressure, air temperature, humidity, and wind speed) agree well with observations with decreasing accuracy in that sequence. This is basically true for the surface radiation fluxes as well, while the turbulent heat fluxes are rather unreliable. On the other hand, some errors such as the turbulent heat fluxes compensate their errors on long-term basis. It was further found that on the long-term some errors tend to compensate (turbulent heat fluxes) and that single models perform significantly better than an ensemble mean depending on the investigated parameter, season, and region.

The availability of RCM model output also fostered the investigation of atmospheric phenomena and processes specific to the Arctic [2227] thereby mostly addressing katabatic winds and related topographical effects (barrier winds, coastal jets) as well as boundary layer features (gravity waves, aerosol distribution). Moreover output from high-resolution climate models can also be directly coupled to glacier mass and energy balance models, as demonstrated for the large ice caps on earth [2831] and is currently underway for Svalbard (Van Pelt et al. [32]).

The aim of this study is to perform climate simulations applying the WRF (Weather Research and Forecasting) model over the Svalbard domain and to validate the skill of the simulations with AWS (Automatic Weather Station) observations at three glaciers on Svalbard. WRF is an advanced mesoscale atmospheric model and is capable to simulate atmospheric processes with high spatial and temporal resolution [33, 34]. Part of the development focussed on optimization of the physical schemes in WRF to improve simulation of the conditions over the Greenland and Antarctic ice sheets. This version is called Polar WRF [35] and is largely based on experience from a precursor model (Polar MM5, [23]). Meanwhile it also accounts for optimum representation of effects over Arctic oceans [36] and Arctic land surfaces [37]. Therefore Polar WRF proved most valuable for investigation of Arctic environments in general [3840] including Svalbard regions [4144].

The latter studies demonstrated the importance of topography on the formation of low level jets at the southern tip of Svalbard (Reeve and Kolsted 2011 [44]), the spatial variability of near-surface meteorological parameters and processes in Svalbard fjords (Kilpeläinen et al. 2011 [41] and Mäkiranta et al. [43]). The latter demonstrated that despite of a high horizontal resolution (1 km) the model was challenged to reproduce, for example, the observed wind direction distribution. This was mainly attributed to insufficient description of local topography and sea ice. The WRF model output performed better for air temperature than for wind as was characterized by RMSE of 1°C and 3 m s−1, respectively.

WRF simulations focussing on Svalbard glaciers are in progress, but validation with relevant AWS measurements has not yet been published. However, a preliminary case study demonstrated the dominance of katabatic winds on glaciers in the Kongsvegen area and their response on local topography and upper air flow [45]. The study showed deficiencies in simulating the boundary layer temperature profile and the jet structure in the fjord which was addressed to insufficient parameterization of surface layer processes.

This paper presents a more comprehensive investigation of the skill of the model. The work aims on identifying improvements towards simulating the atmospheric conditions over glaciers on Svalbard in order to apply the model output for investigation of the mass and energy balance of Svalbard glaciers in a climate perspective. Also we want to find out what horizontal resolution is needed to resolve the local topography. Sensitivity analysis is performed by changing (a) the vertical resolution of the model in the lowest 1500 m, (b) the horizontal resolution to 900 m locally around one of the AWS sites and (c) including varying snow albedo and by changing the physics schemes (turbulence and microphysics). WRF is run with the resolutions 24 km, 8 km for 2 years, in 2.7 km for 3 winter and summer months, and in 900 m resolution and a higher vertical resolution in the lowest 1500 m for July and January. It is forced by ERA-Interim data [46] at the lateral boundary, as well as by sea surface temperature (SST) and ice cover. The model output is validated with data from AWSs placed on three Svalbard glaciers (Kongsvegen, Nordenskiöldbreen, and Vestfonna) measuring air temperature, wind speed, specific humidity, incoming and outgoing short-wave radiation, and long-wave radiation during the period 2008–2010. Precipitation analysis is postponed to later investigations.

2. Model and Measurements

2.1. WRF
2.1.1. Model Setup

In this investigation we use the atmospheric model WRF version 3.2.1 over the Svalbard region [33, 34]. The WRF model is an open-source model (http://www.mmm.ucar.edu/wrf/users/) widely used in research and is useful in ranges 100 of kilometres down to individual cloud scales of 1–5 km. We use the Advanced Research WRF (ARW) dynamical solver, which is a nonhydrostatic model with fully compressible Euler equations and integrates with a third-order Runge-Kutta scheme. Gravity waves and acoustic waves are resolved with a smaller time step. The top of the model is here set at a constant pressure (50 hPa) and with 28 vertical terrain—following Eta—coordinate surfaces in the main simulations. Temperature, horizontal wind, and humidity are calculated between those levels (Arakawa-C grid), adding up to 27 levels. In the main simulations there are 8 levels below ca 1500 m and the lowest two levels at 26 and 90 m. In the sensitivity study with high resolution below 1500 m there are 12 levels below 1500 m, 3 below 100 m, and the lowest two at 11 and 38 m, in all 31 levels.

We use the WRF model in three domains. When simulating the parent domain with 24 km resolution, the initial and boundary forcing from the ERA-Interim archive fields are used every 6 h. The archive represents a reanalysis of the atmospheric conditions every 6 hours, built from ECMWF’s Integrated Forecast Systems (IFS) and available for researchers of member states at http://ecmwf.int/services/archive/. Data of temperature, wind, and relative humidity were used from different pressure levels from 50 to 1000 hPa. Lower initial boundary conditions include soil temperature and moisture at two levels, SST, sea-ice, surface and sea-level pressure, and near-surface wind, temperature, and humidity. The fact that sea ice may cover a fraction of a grid cell was accounted for, which means that the meteorological parameters calculated over water and ice were weighted for each grid according to sea-ice fraction.

Because of numerical instability problems, probably caused by land/sea masking incompatibilities between ERA-Interim and the geographical data from the US Geological Survey (USGS), time evolution of SST and sea ice was handled by doing 10-day consecutive runs. In each run the sea ice and SST were constant. The lower boundary of the soil temperature from ERA-Interim (below 1 m depth) was found too high (of order 20°C), but this is only used at the initialization of each run. In the following time integrations the soil model is fully coupled to the atmosphere part and thus soil temperature dynamically adjusts to the surface temperature forcing. However, the soil temperature is not expected to have a significant influence on 2 m temperature due to the isolating effect of the snow. Results from initial time steps were not used in the analysis. This can be interpreted as allowing the model to spin up for 3 h. The use of ERA-Interim data instead of the NCEP Climate Forecast System Reanalysis (CFSR) data was motivated by the 4D assimilation used in ERA-Interim versus 3D assimilation used in CFSR and the fact that CFSR ends in January 2010.

The largest domain with 24 km horizontal resolution covers the whole Svalbard archipelago. We use a one-way nesting technique meaning that only information from the parent domain is passed on to the inner nests. This approach is common and usually avoids stability problems that might appear when using a two-way nesting procedure [47]. The second nested domain has a grid resolution of 8 km. The innermost domain in the main run has the resolution of 8/3 km (hereafter referred to as 2.7 km). In a sensitivity analysis from one winter and one summer month an inner 8/9 km resolution domain centred at Nordenskiöldbreen was used (hereafter referred to as 900 m). The site was chosen because of its more complex terrain, compared to the other AWS sites. The time step in the model is adaptive and has about 432 seconds for the 24 km nest, about 30 s for the 8 km nest. The 2.7 km nest simulations are performed by downscaling the 8 km simulation afterwards. This keeps the time step as long as 48 s. The 900 s resolution kept the time step at 15 s. Output fields are stored every 3 h. Model domain set-ups are shown in Figure 1 which also shows the improved terrain representation with higher resolution.

The main runs were performed with polar-optimized physics according to the Polar WRF physics [35]. The turbulence and microphysics schemes were interchanged in a sensitivity study to find possible improvements. Note that the polar-optimized WRF should not be confused by Polar WRF itself, which also accounts for varying variable sea-ice thickness and snow thickness over sea ice as well as for seasonally varying ice and snow albedo during the melt period. Comparison between the Polar WRF and the WRF version used in this investigation is envisaged in future studies.

2.1.2. Physical Schemes

The physics schemes are in the main run chosen to reflect those used in the development of Polar WRF [31]. Short-wave radiation is calculated with the Goddard short-wave scheme [48] with 11 spectral bands that account for diffuse and direct solar radiation. Both scattered and reflected components are accounted for. Long-wave radiation uses the Rapid Radiative Transfer Model (RRTM) [49] that is developed to better represent clear-sky conditions. The Grell-Devenyi ensemble [50] convective parameterization scheme is utilized as cumulus scheme. The scheme is of mass flux type and it is run over each grid using different static and dynamic controls. This process gives a total of 144 different convective schemes that are averaged to give feedback to the model. For microphysics the Morrison double-moment scheme [51] is used. The scheme includes six hydrometers: vapour, cloud droplets, cloud ice, rain, snow, and graupel/hail. When predicting the number concentration and mixing ratio (double moment), the particle size distribution is well resolved. This facilitates the calculation of the microphysical processes and cloud distribution. The particle size distribution is treated using gamma functions.

Boundary layer turbulence is parameterized using the Mellor-Yamada-Janjic (MYJ2.5) planetary boundary scheme [52, 53] which must be run together with the Eta surface layer scheme [52, 53]. The MYJ scheme contains a nonsingular implementation of the Mellor-Yamada Level 2.5 turbulence closure model [54] that is active in both the PBL and in the free atmosphere. The Eta surface layer scheme is built on similarity theory [55].

In all WRF set-ups the Noah land-surface soil temperature and humidity model (LSM) [56] is used. The scheme provides sensible and latent heat fluxes to the boundary layer and is able to predict snow cover. The model has 4 soil layers at the constant depths 10, 30, 60, 100 cm. These are converted into snow, with associated snow properties, in deep-snow-covered areas like the glaciers at Svalbard. In the scheme the snow density is set to 0.2 kg m−3, snow albedo to 0.8, and glacier ice albedo to 0.7. On Svalbard the density of snow is higher, between 0.3 and 0.4 kg m−3, because of drifting snow and settling processes [32, 57]. In a sensitivity analysis in the present study monthly climatologic values of the albedo were applied.

In the sensitivity study the WRF single-moment 6-class [58] is utilized. This scheme has the same hydrometeors as the Morrison scheme, but only includes mixing ratios and not particle number concentration. This scheme has been used by Wilson et al. [40] for the Arctic System Reanalysis. For PBL turbulence and surface layer scheme the QNSE or MYNN2.5 scheme is used, respectively. The QNSE scheme (Quasi-Normal Scale Elimination PBL/surface layer) [59] is a    model but uses a new theory for stably stratified regions. It accounts for anisotropic turbulence and internal gravity waves and is not based on the shortcomings of Reynolds stress models in very stable stratifications. This scheme was successfully tested in Kilpeläinen et al. [60] over two fjords in Svalbard. The MYNN level 2.5 scheme (Nakanishi and Niino [61]) includes better interaction of the turbulence with microphysics and radiation in the PBL, compared to the MYJ2.5 scheme.

2.2. Measurements

The model is verified with surface observations from three Automatic Weather Stations (AWSs) located on three different glaciers, namely, Nordenskiöldbreen (NB, 78°41′39′′N, 17°09′22′′E, 530 m a.s.l.), Vestfonna (VF, 79°56′03′′N, 19°11′08′′E, 305 m a.s.l.), and Kongsvegen (KV, 78°46′49′′N, 13°9′22′′E, 537 m a.s.l.). See Figure 1(c) for location of these glaciers, which are representative for the central, northern, and western parts of Svalbard. Close-up maps are shown in Figure 2.

Nordenskiöldbreen (Figure 2(a)) is one of the major outlet glaciers of the ca 600 km2 large Lomonosovfonna ice field draining towards the Billefjorden in the inner part of Isfjorden. The glacier is about 5 km in width and 17 km long and span 193 km2. The AWS is situated in the central flow line of the glacier, which is confined by steep slopes to the north (De Geerfjellet, 1023 m a.s.l. and Flemingfjellet, 1124 m a.s.l.) and to the south (Terrierfjellet, 1211 m a.s.l.), respectively. The direction of the downward slope of the glacier at the AWS site is to southwest. Local topography in the direct vicinity of the AWS is fairly smooth and homogeneous, while on a larger scale the surface is relatively rough due to seracs and open crevasses. Den Ouden et al. [63] and Van Pelt et al. [32] give further details about the topography and the position of the AWS, which is representative for the higher ablation area of the glacier. Descriptions of the site and real-time data are presented on http://www.projects.science.uu.nl/iceclimate/aws/.

Vestfonna (Figure 2(c)) is the second largest ice cap on Svalbard (2400 km2) covering the western part of Nordaustlandet. The AWS is located at the western slope of the ice cap. Local topography in the direct vicinity of the AWS is smooth and homogeneous. The site is characterized by north westerly aspect and located in the ablation area of the ice cap. More detailed topographic features and about the instrumentation may be retrieved from [64, 65].

Kongsvegen (Figure 2(b)) is one of the Svalbard benchmark glaciers where long-term mass balance measurements are performed. The glacier originates at a saddle connecting to the south flowing Sveabreen (800 m a.s.l.) and drains into Kongsfjorden. It covers an area of 102 km2, has a length of 26 km, and slopes have a north westerly aspect. The AWS on Kongsvegen is representative for the average equilibrium line altitude of the glacier, where the surface is virtually smooth and homogeneous. The closest obstructions of the horizon are Dronningfjella (1263 m a.s.l.) about 6 km to the north and Gjerstadfjellet (1006 m a.s.l.) about 5 km to the south. Further details about the topographic features of the glacier and the measurement site are documented (e.g., in [66]).

The surface conditions at all the three sites are characterized by snow during winter and bare ice during summer except on Vestfonna and Kongsvegen where the snow did not melt in 2008. On Vestfonna the ice surface was snow-free in August-September in 2007. At Kongsvegen and Nordenskiöldbreen, snow-ice transitions usually occur in May and September. Surface roughness is generally higher during summer due to differential ice melt and the formation of melt water channels.

The study focuses on validation of the following model output parameters with AWS measurements of air temperature, wind speed, specific humidity, incoming short-wave, and long-wave radiation. Most meteorological variables were measured by similar instruments at the three sites. Kipp and Zonen CNR1 instruments were used to measure incoming and outgoing short-wave and long-wave radiation at all sites. Wind speed and direction were measured by Young Wind Monitor RM 05103. Humidity and air temperature were measured using Vaisala HMP45A probes mounted within naturally ventilated shields (Vestfonna and Nordenskiöldbreen) while an artificially ventilated device (Rotronic Hygroclip) was used at the Kongsvegen site. The used humidity sensors are calibrated to measure relative humidity with respect to water and the data were converted for reference to ice based on the ratio of saturated vapour pressures. Mean hourly data were stored onto Campbell Sci. CR10X and CR23X loggers and the sites were usually visited in spring for maintenance.

Meteorological data collected by the AWS at Nordenskiöldbreen used in this study are from March 2009 to March 2010. Unfortunately, the temperature readings were erroneous during the summer of 2009 (July to November). This also affected the humidity, long-wave radiation, and snow height readings. For this period the logger temperature was used. To that end the logger temperature was correlated to the air temperature using a linear regression (). Using this method, the estimated air temperatures during periods with high solar input may be overestimated by up to 5°C. All the measurements at Nordenskiöldbreen are observed approximately 4.5 m above the surface. Observations at Kongsvegen and Vestfonna are from January 2007 to December 2010 and May 2007 to April 2009, respectively. On Vestfonna all of the parameters are measured at the height 1.8 m except for the wind speed that is measured at the height 2 m. At Kongsvegen air temperature and humidity are measured about 1 m below the wind measurements corresponding to a height of 3.5 m above the glacier ice. Note, however, that in glacier environments the height of instruments above the surface inevitably varies throughout the year due to accumulation and melt of snow. To exemplify this for the Kongsvegen site, air temperature was measured about 0.9 m above the snow (May 2008) and 3.0 m above the ice surface in September. The initial snow depth at Nordenskiöldbreen in March 2009 was 0.85 m. The varying measurement heights were not corrected for within this study because of the lack of knowledge of the near-surface temperature profile. However, one can qualitatively argue that the 2 m temperature is higher than the recorded temperature when the measuring height is less than 2 m, as a surface inversion is common over the snow and ice surfaces.

Measurements in the harsh Arctic environments impose some further constraints and enhanced treatment of the data. At Kongsvegen, for instance, values of temperature outside the range −40°C to +20°C and humidity data (visual outliers) were removed. High temperature peaks during the summer at Nordenskiöldbreen, thought to be associated with errors caused by short-wave radiation warming, were visually identified and the upper limit of +5°C was used to sort out these occasions. Wind speeds less than 0.5 m s−1 were removed due to the inability of the instruments to correctly measure these low winds and also because of the possibility that the equipment might have been covered in ice. Short-wave radiation data in excess of 1000 W m−2 were removed and potential snow on the up-facing radiometers was sorted out by values albedo higher than 0.95. Long-wave radiation data at Kongsvegen were controlled by comparison with values calculated according to König-Langlo and Augstein [67] and the 0°C threshold for calculated surface temperature (melting ice), respectively. Remaining gaps in the record were closed by cross correlation using undisturbed data measured concurrently at Ny-Ålesund and at the glacier.

All mentioned data were provided by different research groups as detailed in the acknowledgments.

3. Methodology

The comparison was made using model data from one out of the four grid points that surrounded the AWS station. The closest grid point may deviate in height from the AWS and therefore not be a suitable for comparison. When choosing grid point corresponding to an AWS the correlation to the observations and the topography was regarded and to some extent also the bias. If there was no significant difference between the correlations among the grid points, the point with the best coincident height was used. It appeared that the wind speed was the most sensitive variable to location. The resulting topographical heights are for (1) Kongsvegen: 278 m a.s.l. (ERA-Interim), 424 m (24 km), 639 m (8 km), and 607 m (2.7 km), (2) Nordenskiöldbreen: 380 m (ERA-Interim), 707 m (24 km), 321 m (8 km) and 648 m (2.7 km), and 607 m (900 m), and (3) Vestfonna: 260 m (ERA-Interim), 188 m (24 km) and 244 m (8 km). The grid points in the 2.7 km resolution consist of 100% of glacier, whereas the 24 km grid points at least are of 53% covered by ice (Nordenskiöldbreen).

In many areas the temperature decreases with height and one can for climate studies use constant lapse rates. The latter is not a good approach when considering the effect of the diurnal cycle and the complex nature of stable boundary layers, which is very common in Arctic areas. Particulary during the winter time when inversions are common, a temperature correction accounting for the altitude differences could even increase the error. In a test of Polar WRF carried out over Greenland by Hines and Bromwich [35] a correction of a decrease in temperature of 0.0071 K m−1 [68] was used. In Pohjola et al. [69], a mean lapse rate of −0.0044 K m−1 was found for the area of study. Here the temperature was both analysed without a temperature correction and with the −0.0044 K m−1 lapse rate. The maximum altitude error between the AWSs and the WRF model or ERA-Interim grid in this investigation is 250 m, which corresponds to approximately 1.8°C. The varying heights of the AWS instruments measuring temperature are not accounted for.

To investigate if the WRF downscaling improved the modelled data compared to ERA-Interim output, the WRF standard output for 2 m temperature, 2 m specific humidity, 10 m wind, and incoming short-wave and long-wave radiation were compared with mean hourly observations from the AWSs every 3 hours. The skills of WRF with polar physics output were evaluated for the four different resolutions: 24 km, 8 km, 2.7 km, and 900 m). The 24 km and 8 km resolution runs were performed for the period April 2008–March 2010, covering 2 years. Because of limited computer resources 2.7 km was only simulated over 3 months in summer (June, July, and August = JJA) 2009 and winter 2009-2010 (December, January, and February = DJF) driven by 3 hourly boundary conditions of the parent 8 km domain. Sensitivity simulations were performed for these summer and winter months. Additional runs with an improved horizontal resolution of 900 m were simulated for July 2009 and January 2011.

In the sensitivity runs with different physics schemes the same physics schemes were used in all nests from 24 km to 2.7 km (and 900 m). Cumulus scheme was turned off for horizontal resolutions below 8 km. In the simulation with higher vertical resolution we expect that the stratification and wind profile of the boundary and surface layer will be simulated more accurately. In a modelling study by Söderberg and Parmhed [70] of katabatic wind jets over an Icelandic sloping glacier with complex surrounding terrain the vertical resolution was set to 2 m in the lowest 21 m while the Mellor Yamada original PBL scheme [54] was used. This was motivated by the large wind gradients in the surface layer. Also stable surface layers are very shallow.

Simulated data were evaluated using correlations coefficient (), mean difference (Bias), mean absolute error (MAE), and skill score coefficient (SSC). Following [71], we defined the latter as the mean of the absolute difference between the modelled data and observations which is normalised with standard deviation of the observations :

4. Results

In the following section the input from ERA-Interim and the output from our runs with WRF 24 and 8 km data during the period April 2008–March 2010 are presented. Because of limited computer resources the fine 2.7 km resolution simulations were only performed during a summer and winter period where data from two out of three AWSs were available. Therefore the periods JJA in 2009 and DJF in 2009-2010 were chosen. This had the consequence that Vestfonna is excluded in that analysis. It was also expected that the finer resolution was more important for the other sites with more complex topography. In the statistical comparisons with 2.7 km the same period is used for 8 km.

4.1. Impact of Improved Horizontal Resolution

ERA-Interim gives overall very good correlations for temperature, probably due to the assimilation of data from synoptic weather stations like Longyearbyen and Barentsburg. The 2-year simulations with the WRF model, compared to ERA-Interim, show only partly improved biases and correlations with the AWS observations. Scores improved hardly with increasing resolution (Table 2). For wind speed it is clearer that the downscaling with WRF increases the ability to simulate the wind speed. ERA-Interim is not able to capture the variability of the wind and the correlation is very low.

The observed temperature at Vestfonna, Kongsvegen, and Nordenskiöldbreen varies in the ranges −35 to 6°C, −34 to 10°C, and −33 to 9°C, respectively. These ranges are underestimated by ERA-Interim but improve with finer resolution, at least to 8 km. The diurnal temperature cycle (not shown) is underestimated in the WRF model and becomes smaller with better resolution. Comparing all three WRF resolutions (winter and summer) there is no significance in that the correlation in temperature increases or that the bias decreases with resolution (Table 3). However, the bias can depend on the actual height of the AWS in relation to the model points. Correcting for this (-corr) by applying the temperature lapse rate −0.0044°C/m (discussed in Section 3) the bias, MAE, and SSC are partly improved. This is more important for the ERA-Interim and coarse resolution simulation data (Table 2 and for summer and winter in Table 3). It appears that the WRF model fails especially when simulated long-wave radiation has large errors (Figure 3(a)). Underestimation of incoming radiation gives large underestimation of the temperature in the model because of the energy balance. Also there is a tendency that the lowest observed temperatures are not captured and that the variation of the summer temperature is underestimated (Figure 3(b)). The humidity is closely coupled to the temperature and follows largely the same reasoning as for temperature (for statistical details confer Table 4). Humidity was not measured during the summer at Nordenskiöldbreen (see Section 2.2).

The wind speed is more problematic to analyse than 2 m temperature. The anemometers at height are all below 10 m at the AWSs and will presumably underestimate the 10 m wind. Therefore one can estimate a 10 m wind speed from the original AWS wind speed based on, for example, Monin-Obukhov theory. Considering all uncertainties of the processes in the stable surface layer we base our analysis on a very simple correction. An increase of the wind speed by 20% at Kongsvegen, 15% at Nordenskiöldbreen, and 25% at Vestfonna was applied. This relatively rough correction is based on a roughness length,   mm, applicable on a snow surface [72] and stable stratification with an Obukhov length of 50 m at Kongsvegen and Nordenskiöldbreen and 100 m at Vestfonna (higher wind speeds). As we do not have any information of the actual atmospheric stratification (apart from the simulated) we start with this approximation. The correction will not change the correlations but the biases, MAEs, and SSCs. The variability of the corrected wind is relatively well captured but the fine resolution has a tendency to overestimate the highest wind speeds. This is considered as an artefact of the applied constant correction as expected for more neutral stability at high wind speeds. The statistical evaluation is presented in Table 2 and for summer and winter in Table 5. It appears that increasing resolution does not significantly improve the scores at Nordenskiöldbreen and Kongsvegen. At Vestfonna with much less complex topography increased resolution has a higher positive impact (Table 1). This suggests that 2.7 km resolution is not sufficient to resolve the important wind speed determining features of the local topography in more complex terrain. The finding here of better results with coarser resolutions might hence be an artefact, being caused by randomness in the location of the grid point and the neighbouring topographic structure in that resolution. There is no significant dependency on certain conditions found when the simulated wind speed is deviating from the observed.

The wind direction is however rather well described. In Figure 4 wind roses are shown for Nordenskiöldbreen and Vestfonna (Kongsvegen data were not available). It is clear that the WRF model is able to simulate katabatic or other slope winds, even though we cannot not say anything about its vertical profile and the quantitative values. Note that at the considered Nordenskiöldbreen grid point is located south of the AWS where downward slope is more to the west. Due to the complex topography also the wind direction is more spread out in the simulation.

Due to the high-latitude locations of the AWSs short-wave radiation is limited to the summer period, when there is also midnight sun. The incoming short-wave radiation, both simulated and observed, seldom exceeds 600 W m−2. The statistics are given in Table 6 for summer and winter and it can be concluded that 8 km resolution improves the correlation. However, the bias and errors are not necessarily improved. The overestimation is largest for the finest resolution which indicates a nonlinear problem. The effect of increasing the complexity of the topography affects the cloudiness and if not sufficiently well described by the resolution it may induce more errors in the cloud cover, compared to a coarser resolution. Here it is evident that the model underestimates the cloud cover. The radiation from the sun also penetrates the clouds. The simulated incoming short-wave radiation amount during cloud covered conditions is also very important for radiation balance [73]. The performance of the WRF model and the polar-optimized short-wave radiation scheme during overcast was however not analysed in this investigation.

For the energy balance at the surface the outgoing short-wave radiation is equally important as the incoming short-wave radiation. The outgoing part is determined by the surface albedo, which in the present setup of the model is set constant to 0.8 for snow surfaces. Figure 5 shows the measured albedo at Kongsvegen, defined as the ratio of outgoing and incoming short-wave radiation, together with the albedo in the 2.7 km resolution WRF. Typically the albedo decreases rapidly during the melt season, due to snow metamorphism, enhanced liquid water content, enrichment of impurities, and dust deposition of for instance black elemental carbon [74, 75]. These effects are implemented in Polar WRF, but not in the WRF version used in the present investigation. In this model version albedo over Kongsvegen is fluctuating due to snow melt thereby episodically exposing the glacier ice surface. Note that according to the measurements new snow has albedo above 0.9, which is not captured by the model. Thus it is clear that varying albedo is not yet realistically represented which is a key issue when simulating ice/snow-air exchange on the other hand. The built-in function in WRF to use monthly climatology values of the albedo was applied and tested in the summer period June–August 2009. However, this study did not yield a significant improvement compared to using the table values used as standard in WRF. At Nordenskiöldbreen it varied between 0.08 and 0.55 and at Kongsvegen between 0.07 and 0.7. The resulting 2-meter temperature was however not significantly affected, only by lowering it by 0.1°C in mean. The diurnal temperature variation was reduced by less than 0.2°C as well. It should here also be stated that the observations might not be completely representative to the model grid.

Incoming long-wave radiation is very dependent on whether the sky is clear or cloud covered. Particularly low clouds increase the incoming long-wave radiation, up to 100 W m−2. Table 7 shows the statistical evaluation during summer and winter. As in the example of short-wave radiation the 8 km resolution gives somewhat better scores than 2.7 km. The former discussion regarding topography not resolved by the model also applies here. The influence from clouds is well shown in Figure 6, where the four nodes define the ability of the model to capture or miss the presence or nonpresence of clouds. There are occasions where the model both over- and underestimates the cloud cover. It appears that the WRF model simulates the average cloud cover quite well but with a negative bias of order 10 W m−2 in the winter and smaller negative or positive bias in the winter.

4.2. Sensitivity Studies

A number of sensitivity simulations were performed for 1- or 3-month summer and/or winter periods. The simulations are summarized in Table 8 together with CPU times. It is evident that the WRF6 microphysics requires less CPU time compared to the more complex Morrison double-moment scheme. The QNSE and MYNN2.5 schemes require more computational resources than the MYJ2.5 scheme. Further it is seen that the finer horizontal resolution simulations are more computer expensive than increasing the number of vertical layers in the lowest 1500 m by 4.

The sensitivity analysis of the different simulations, except the albedo analysis, are presented in the next two subsections.

4.2.1. Impact of Parameterisation Scheme

The statistical evaluation of the different physics schemes is given in Tables 9, 10, 11, and 12. For temperature (Table 9) it is evident that the correlation coefficient is not remarkably changed, but the QNSE scheme is slightly better in winter, compared to the polar physics. MYNN2.5 is better correlated to the observations in summer and simulation of air temperature is rather insensitive to the choice of microphysics scheme. The bias and errors show some improvements and follow about the same patterns with the preference to QNSE in winter and MYNN2.5 in the summer. In Figure 7 it is shown that the different schemes do not change the distribution to a large extent. Both low and high temperatures are often underestimated, but QNSE scheme provides the best distribution. The mean diurnal temperature variation is underestimated in all simulations (not shown), except for the MYNN2.5 scheme in summer at Nordenskiöldbreen. Overall, the MYNN2.5 scheme provides the largest and best diurnal variation.

Also for the wind speed, the choice of microphysics is of minor importance. The QNSE scheme improves the correlation in winter at Kongsvegen but worsens it during the summer. The correlations are very low, between 0.27 and 0.55. The MYNN2.5 scheme is either worsening or only slightly improving the correlation to observations. The errors and biases should not be analysed in too much detail because of the height corrections made. But although the observed wind is corrected upwards the model still overestimates the mean wind. It is often an improvement to use the QNSE or MYNN2.5 schemes instead of the MYJ2.5 scheme, MYNN2.5 giving the lowest biases. For the wind distribution in Figure 8 (omitting all values below 0.5 m s−1) it appears that the simulated distributions are too flat, that is, the model produces too many cases with low and high winds. There is no clear scheme that outperforms the others but the MYNN2.5 scheme yields less frequent cases with strong winds and is in average closer to the observations.

Long-wave radiation correlation is worsened with the microphysics scheme except in summer at Kongsvegen and biases and errors follow about same pattern. The turbulence schemes, however, improve in mean correlation, biases, and errors with the QNSE scheme being slightly superior to the MYNN2.5 scheme, especially in winter.

The short-wave radiation is, as expected, most affected by the microphysics due to the processes that determine the cloud extent and thickness and how much of the moisture leaves the clouds via precipitation. The correlation is increased by the WRF6 scheme (Table 11). However, the bias increases and errors are not improved. The PBL schemes somewhat improve the biases, probably because of changed humidity fluxes over the open sea.

It appears that the use of the simpler microphysics scheme WRF6 in average reduces the skill of the WRF model on average, compared to the Morrison scheme. The improved correlation for short-wave radiation was cancelled by the worsened biases and errors and the worsened skills for long-wave radiation. The shorter CPU times with the scheme might not be worth in relation to model performance. Regarding the turbulence schemes both QNSE and MYNN2.5 improve the model skills, compared to the MYJ scheme. However, finding the best of them may require further investigations. The tendency here is that QNSE works better in winter and MYNN2.5 in the summer.

4.2.2. Impact of Improved Horizontal and Vertical Resolution

Because of limited computer resources the high-resolution simulations were only performed for July 2009 and January 2010 and the 900 m resolution only over Nordenskiöldbreen. Also the higher vertical resolution caused more numerical instability (when trying to simulate the other winter and summer months) than the other simulations, although an adaptive time step was used (the QNSE turbulence scheme instead of the MYJ2.5 scheme improved the stability but was not analysed). The higher vertical resolution in the boundary layer improves the correlation of the simulations to observations of some of the meteorological parameters, temperature (Table 13), wind (Table 14), and incoming long-wave radiation (Table 15). The 900 m resolution increases the correlation only for long-wave radiation. For temperature and wind the bias and MAEs are not significantly improved (or worsened). Time series of the high-resolution simulations only occasionally changes the temperature (not shown). In Figure 9 the distribution of temperature at Kongsvegen and Nordenskiöldbreen is shown. The simulated temperatures are bias-corrected based on the January and July bias, respectively. The observations have a double-peak in the winter and a single one in the summer, while the simulated distribution has double peaks in both seasons. The distribution is only slightly improved by the increased vertical resolution. The 900 m resolution shows more significant improvement. The mean diurnal variation is underestimated by the simulations except for Nordenskiöldbreen in winter. The higher vertical resolution decreases the variation, whereas the 900 m resolution increases the variation. Our estimate of the observed 10 m wind speed is rough and we cannot expect that the bias and errors would be improved by the higher vertical resolution. The 900 m resolution however lowers the simulated wind and lowers the errors. The wind distribution (Figure 10) is not much improved by the vertical resolution but the horizontal resolution better simulates all wind speeds above 7 m s−1. The wind direction is very dependent on horizontal resolution as it makes it possible to find grid points with the correct terrain slope direction. In Figure 11 and comparing with the observations and with 2.7 km resolution in Figure 4 it is evident that the wind direction is very much improved even if the side winds are still slightly exaggerated compared to the observations. Even when only increasing the vertical resolution, the wind direction becomes less spread compared to the main runs but the principal direction is too easterly.

The improvement of the long-wave radiation in correlation, bias, and MAEs for higher vertical resolution can probably be traced back to better resolved cloud formation in the boundary layer (Table 15). Increased horizontal resolution does not improve the correlation in winter but otherwise the bias and errors are improved with this configuration. In Figure 12 time series in January and July of incoming long-wave radiation at Nordenskiöldbreen are shown. The simulated long-radiation relatively well follows the observed fluctuations but underestimates the radiation in periods with less than 50 W m−2. These low values are associated with clear sky or missed clouds in the model. Also when the radiation is high in both observations and in simulations the long-wave radiation is somewhat underestimated, except in the end of July. Most of the time the difference is very small between the simulations, but occasionally either the higher vertical or the higher horizontal resolution captures the changes in incoming long-wave radiation. For instance, the increased vertical resolution captures the clear sky in the beginning of January and the low clouds on 12-13 July, whereas the 900 m resolution captures the clouds on 17-18 July. Incoming short-wave radiation (Table 16) is improved with the higher vertical resolution but not as significantly as the long-wave radiation. The correlation for the 900 m simulation decreases but the bias and errors are improved.

To summarize, the higher vertical resolution improves all the simulated parameters, except the distribution of temperature, wind speed, and direction. Regarding the latter the 900 m horizontal resolution performs best. The changes for the other parameters are though rather small and in view of CPU time the higher vertical resolution seems to be the most economic change of the grid resolution. The 10 m wind direction is slightly improved (less spread) also with the increased vertical resolution. This suggests to potentially use a grid point with a terrain slope reflecting the actual observation site. However, when testing surrounding grid points in the 2.7 km resolution none of them gave better wind directions compared to the grid point used in the investigation, based on the correlation to the wind speed. Further, we have used a rather modest increase of vertical resolution of the PBL. It is possible that there would be even better scores of the meteorological surface layer variables with finer vertical resolution. As mentioned in Section 3, Söderberg and Parmhed [70] used 2 m vertical resolution in the lowest 10 s of meters.

5. Discussion and Conclusions

Results are presented from simulations performed over Svalbard by the mesoscale climate model WRF 3.2.1 for 2 years from spring 2008 to spring 2010. The ERA-Interim reanalysis was dynamically downscaled with the WRF model using the physics setup used in Polar WRF. The polar-optimized WRF model simulated the climate in 3 domains with the resolutions 24 km, 8 km, and 2.7 km. The inner nest was only simulated for summer 2009 and winter 2009-2010. Also sensitivity simulations were performed with one alternative microphysics scheme, WRF6 (used in the Arctic Re-Analysis), and two PBL/surface layer schemes, QNSE (developed for very stable conditions) and MYNN2.5 (with turbulence interaction with clouds). Additional simulations were performed for July 2009 and January 2010 considering higher vertical resolution in the lowest 1.5 km and 900 m horizontal resolution, respectively. The investigation focuses on validating the following surface variables: temperature, wind speed, specific humidity, and incoming short-wave and long-wave radiation. The WRF data are compared with observational data received from three automatic weather stations on Svalbard, Kongsvegen, Vestfonna, and Nordenskiöldbreen.

Downscaling the ERA-Interim analysis with the WRF model improved the simulated wind speed and to some extent the temperature. The correlation of the ERA-Interim temperature with the observations was from the beginning very high, being around 0.9, probably due to the assimilation of data from nearby SYNOP stations. Temperature changes are thus less sensitive to the topography than the wind speed.

Increased resolution did not always give better results, probably due to important unresolved topographic and heterogenic effects even for the 2.7 km resolution. The correlation of the wind speed in the 2.7 km resolution ranges between 0.32 in winter (Kongsvegen) and 0.55 in summer (Nordensköldbreen). The terrain is rather complex at both glaciers. It was examined if the lack of improvement from increasing the resolution from 8 km to 2.7 km may also be attributed to possible shortcomings with the downscaling method. The downscaling to the 24 km and 8 km grid is performed every time step (c. 7 min), meaning that every calculated field in the outer domain is used as boundary condition to the inner domain. The downscaling from 8 km to 2.7 km, done to reduce the CPU time, only provides boundary conditions each 3 h. Thus the 2.7 km nest lacks the transitional information between the 3 h inputs. However, a sensitivity study performed in the winter period showed however that this was of minor importance. The correlations for wind and temperature and the bias for wind at Vestfonna, with smoother topography, were higher than the other observations and improved with increased resolution. The radiation is very much influenced by clouds and at least the low-level clouds appear to be underestimated, leading to periods of underestimated incoming long-wave radiation and overestimated incoming shortwave radiation. As a consequence the winter temperature was underestimated (when there is no short-wave radiation). During the summer months the bias was smaller or even positive due to the erroneous simulated shortwave radiation.

Higher vertical resolution in the boundary layer improved the simulations of radiation and wind speed, probably due to more correct vertical profiles of wind speed and moisture transport from the open sea. Temperature was not significantly improved. Higher horizontal resolution over Nordenskiöldbreen did surprisingly not improve the wind speeds, but the direction and also frequency distribution of temperature and wind. It requires more computer resources to increase the horizontal resolution in relation to the vertical. The room for improvement regarding resolution would thus favour the vertical. The vertical resolution may very well be increased more than 8 layers below 1500 m and would probably increase the scores.

The sensitivity study of the physics schemes revealed that for the Kongsvegen and Nordenskiöldbreen sites the Morrison microphysics scheme is recommended as is used in the Polar WRF. It is also more sophisticated than the WRF6 scheme giving improved incoming long-wave radiation, due to better cloud simulation. Both the QNSE and the MYNN2.5 schemes showed improvements compared to the Polar WRF scheme MYJ2.5. The QNSE scheme was slightly better than MYNN2.5 in winter and the opposite in the summer. When performing the high vertical resolution simulations it showed out that the model was suffering from numerical instability dependent on the PBL scheme. It might thus be crucial to use a certain PBL scheme when increasing the vertical resolution.

Overall, this investigation proves a good ability of the WRF model to simulate present climate on Svalbard glaciers, except for wind speed where correlations are as low as . However, the importance of having physics schemes suitable for polar environments and fine vertical and horizontal resolution cannot be overemphasized. The study points out that simulations mainly could benefit from better descriptions of surface roughness length, and seasonal and shorter-term variation of glacier albedo (snow and ice), which is partly realized in Polar WRF already. Short fluctuations in the surface albedo are associated with new snow and high albedo, melting snow with lower albedo while even older snow or bare ice are strongly affected by dust deposition darkening the surface. In the measurements the albedo ranges between 0.2 and 0.95, whereas in Polar WRF it is varying between 0.5 and 0.82. Even though the 2 m temperature was proved to not be very sensitive to albedo it is of major importance for the glacier/snow surface energy balance.

Main conclusions regarding the WRF model performance based on observations at Kongsvegen and Vestfonna are as follows. (1)For future simulations with the WRF model over Svalbard glaciers, (a)polar WRF physics should be used, except regarding the PBL scheme, where QNSE or MYNN2.5 better describes distribution of temperature and wind speed, the QNSE scheme being better in winter; (b)8 km resolution is not sufficient to model fine scale variations of the field but may be sufficient for Vestfonna due to a much simpler topography; (c)finer vertical resolution in the lower part improves the model score, but probably there is need for even finer resolution than that used in this study (3 levels below 100 m and lowest at 11 m). (2)Improvements are needed regarding albedo, turbulence, and the related surface roughness.

Future possible investigations include precipitation validation and quantitative wind direction evaluation. The influence of terrain gives realistic enhancements of the precipitation in relation to the surrounding sea, but it remains to quantitatively evaluate the model with observations and whether the Morrison scheme outperforms the WRF6 scheme also for precipitation amount. Evaluating simulations with even finer vertical resolution and with the QNSE and MYNN2.5 scheme is also very interesting. Regarding observations there is firstly a need for more and consistent long-term meteorological observations at these glaciers. Secondly, alternative parameters may be considered. Thus turbulence measurements over glaciers would help to find out if the model performs well with respect to the turbulent flux components of the energy balance. Momentum flux measurements can also give information of surface roughness.


The authors thank their research groups and funding agencies for providing the observation data and related meta information. Thus, the Vestfonna data were acquired within IPY-KINNVIKA project expeditions funded by the Swedish Research Council. The data were provided by Regine Hock (Geophysical Institute, University of Alaska), Rickard Pettersson, and Ulf Jonsell (Swedish Polar Research Secretariat, Stockholm, Sweden). Thanks are due to Ulf Jonsell also for comments on the paper. The Vestfonna data are also available through http://www.smhi.se/ecds/Search-data. The achievement of data from Nordenskiöldbreen station was financed by The Netherlands Organisation for Scientific Research (NWO) and the Swedish Science Council as part of the IPY-GLACIODYN project. The glacier field work at Vestfonna and Nordenskiöldbreen was partly supported by Norwegian Polar Institute (NPI). Kongsvegen data were collected within field work funded by ARCFAC V 2007 (EC Contract no. 026129-75). Fieldwork on Kongsvegen was partly supported by NPI (J. Kohler, Tromsoe). Kongsvegen data analysis was supported by the Austrian Science Fund (FWF, Grant I 369-B17) and ESF, that is, the SvalGlac project (http://svalglac.eu/). NPI is acknowledged for the DEM data. The WRF simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX). ECMWF ERA-Interim data used in this study have been obtained from the ECMWF data server. Work on the paper itself was a joint effort by ESF-SvalGlac members and thereby was also funded by the Nordic council of ministers (SVALI publication no. 6). The authors finally acknowledge two anonymous reviewers for their input to improve the original paper.