The surface energy balance at the Svalbard Archipelago has been simulated at high resolution with the Weather Research and Forecasting Model and compared with measurements of the individual energy fluxes from a tundra site near Ny-Ålesund (located north of Norway), as well as other near-surface measurements across the region. For surface air temperature, a good agreement between model and observations was found at all locations. High correlations were also found for daily averaged surface energy fluxes within the different seasons at the main site. The four radiation components showed correlations above 0.5 in all seasons (mostly above 0.9), whereas correlations between 0.3 and 0.8 were found for the sensible and latent heat fluxes. Underestimation of cloud cover and cloud optical thickness led to seasonal biases in incoming shortwave and longwave radiation of up to 30%. During summer, this was mainly a result of distinct days on which the model erroneously simulated cloud-free conditions, whereas the incoming radiation biases appeared to be more related to underestimation of cloud optical thickness during winter. The model overestimated both sensible and latent heat fluxes in most seasons. The model also initially overestimated the average Bowen ratio during summer by a factor of 6, but this bias was greatly reduced with two physically based model modifications that are related to frozen-ground hydrology. The seasonally averaged ground/snow heat flux was mostly in agreement with observations but showed too little short-time variability in the presence of thick snow. Overall, the model reproduced average temperatures well but overestimated diurnal cycles and showed considerable biases in the individual energy fluxes on seasonal and shorter time scales.
The Arctic region has experienced a larger warming than the global average over the last decades (ACIA 2005; Chylek et al. 2009) and is also the region for which the largest future warming is expected (Stoker et al. 2013). Among the mechanisms contributing to this Arctic amplification is the reduction in surface albedo with melting of snow and ice (Serreze and Barry 2011; Taylor et al. 2013), which is closely related to the surface energy balance (SEB). The thawing of permafrost, which has the potential to release large amounts of carbon to the atmosphere (McGuire et al. 2009), is another important atmosphere–cryosphere interaction influencing the global climate system. Our ability to realistically simulate these processes is limited by our understanding of the SEB and the quality of our models, both of which can be improved by a combination of modeling and observational efforts. Such observations are, however, currently scarce in the Arctic, and there is often a scale gap between the horizontal extent for which the observations are representative and the resolution of the global models. High-resolution, regional modeling is therefore necessary to evaluate the ability of current models to capture the physical processes governing the SEB.
Situated in the high Arctic and consisting of a combination of tundra with underlying permafrost and glaciers, the Svalbard Archipelago is well suited for studying the Arctic SEB. A large number of observations are available within this region, including measurements of all components of the SEB from a permafrost site located north of Norway near Ny-Ålesund (Westermann et al. 2009), but with its large horizontal inhomogeneity the above mentioned scale gap is a challenge.
In this study, we use the Weather Research and Forecasting (WRF) Model to simulate the SEB on Svalbard for one whole year. This model has been used in previous studies of surface-layer and boundary layer processes in this region (Kilpeläinen et al. 2011, 2012; Makiranta et al. 2011; Claremar et al. 2012; Mayer et al. 2012). In addition, a polar-optimized version of WRF (Polar WRF) has been developed and tested for Arctic conditions (Hines and Bromwich 2008; Bromwich et al. 2009; Hines et al. 2011) and is used for a high-resolution Arctic-system reanalysis (Bromwich et al. 2012). The aim of our study is to use detailed in situ measurements from Svalbard (Westermann et al. 2009; Schuler et al. 2014) to 1) assess the ability of this model to capture the characteristics of the SEB in this region, especially in relation to permafrost and glacier modeling, and 2) to identify model weaknesses for which improvements are needed to capture the important factors controlling the SEB.
2. Model and measurements
a. Model description
The WRF Model, version 3.4.1 (Skamarock et al. 2008), is used in this study, with three one-way nested domains at 9-, 3-, and 1-km horizontal grid spacing and no less than 20 grid cells between each boundary, including the relaxation zone of 4 cells. The domains approximately follow earlier studies of the Svalbard region by Kilpeläinen et al. (2011, 2012). In the vertical direction, 35 terrain-following layers are used, with the lowest layer centered at approximately 27 m above the ground and a total of 7 layers in the lowermost 1000 m.
For the land surface model (LSM) in WRF we use the Noah LSM (Chen and Dudhia 2001) with four soil layers centered at 0.05, 0.25, 0.7, and 1.5 m below the surface. For cloud microphysics we use the Morrison two-moment scheme (Morrison et al. 2005). The longwave radiation and shortwave radiation are calculated using the scheme from the Rapid Radiative Transfer Model for GCM Applications (RRTMG; Iacono et al. 2008) and the (old) Goddard scheme (Chou and Suarez 1999), respectively, except for the innermost domain, for which the RRTMG scheme is used for both longwave and shortwave radiation.1 The planetary boundary layer (PBL) is simulated with the Mellor–Yamada–Janjić scheme in conjunction with the Eta surface-layer scheme (Janjić 2002). Subgrid cumulus clouds are parameterized only in the outermost domain using the Grell–Devenyi ensemble scheme (Grell and Devenyi 2002). The physical options largely follow the ones used in Polar WRF for Arctic land conditions (Hines et al. 2011), even though we do not use the fully modified Polar WRF Model here.
Certain parts of the model with special relevance for the SEB will be described in somewhat more detail. These are the simulation of sensible heat flux Qh and latent heat flux Qe (hereinafter referred to as the turbulent fluxes), surface albedo, ground heat flux Qg, and the simulated residual (Res).
The turbulent fluxes are calculated in the Noah LSM with bulk transfer formulas (Chen et al. 1997), using a common exchange coefficient Ch for heat and moisture, calculated with Monin–Obukhov similarity theory in the Eta surface-layer scheme. On the basis of the formulation of Zilitinkevich (1995), the roughness length Z0t for heat and moisture is related to the roughness length for momentum Z0m:
Here k is the von Kármán constant (=0.4), Re is the roughness Reynolds number, and Czil is an empirical coefficient (=0.1 in our simulations). For the tundra site, Z0m is between 0.1 m (snow free) and 0.003 (snow covered).
The surface albedo is calculated by weighting the snow albedo and the background (snow free) albedo by the snow-cover fraction. The snow albedo is in turn calculated by combining a fixed land-cover-dependent value with a snow-albedo scheme that accounts for the age of snow (Livneh et al. 2010). The latter uses a value of 0.85 for fresh snow and is gradually reduced with time since last snowfall. In our case, this yields maximum snow albedos of 0.835 and 0.75 for glacier and tundra surfaces, respectively.
Flux Qg is calculated from the temperature gradient between the surface (skin temperature) and at the center of the first soil layer. With Noah treating snow as a part of the first layer, a combined conductivity for the snowpack and the soil is used in the presence of snow (Ek et al. 2003).
The residual of the fluxes in the simulation mainly arises from the following processes, which are accounted for by the model but are not given as standard output: the latent heat release when rain freezes at the surface, energy released/consumed as a result of different temperature of precipitation and the surface, and energy consumed by snowmelt. In addition, a (small) imbalance in the calculated SEB exists in the Noah scheme.
b. Simulation strategy
As lateral boundary conditions, we use the ERA-Interim reanalysis, downloaded with 0.5° horizontal resolution and 6-h time steps (Dee et al. 2011). Because initial test simulations of summer 2008 revealed too-high sea surface temperatures (SST) along the west coast of Svalbard in this reanalysis, the SST and sea ice fraction were replaced with daily fields from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA; Donlon et al. 2012), which, from February of 2009, are also used in ERA-Interim. Land-cover data from the U.S. Geological Survey in the model were replaced with those of Nuth et al. (2013), because the former characterized too-large areas as ice covered, and the surface elevation data were replaced with data from the Norwegian Polar Institute.
To spin up soil temperatures, the two coarser domains were initiated with no snow and 0.5°C temperature in all soil layers at 1 September 2007, and data from the middle domain (3 km) were used to initiate the inner domain (1 km). After another 2-day spinup for the atmosphere, the full model was run from 15 March 2008 to 15 March 2009 without nudging or reinitialization. This was done to avoid validating the model against observations that are already assimilated into the reanalysis, assuming that the full domain (about 1000 km by 1000 km) would be small enough to not drift substantially from the reanalysis controlling the boundaries.
A set of measurements from across Svalbard is used to validate the model (Fig. 1). The main validation sites are Ny-Ålesund (red dot), where the full SEB has been measured together with a large number of ancillary measurements, and Austfonna, Norway (yellow dot), where all four components of the radiation budget together with meteorological observations are available. In addition, five other locations with synoptic measurements are used for model validation (white dots).
The Bayelva climate and soil monitoring station is located about 2 km west of the village of Ny-Ålesund (78°55′N, 11°56′E) at about 20-m elevation. It is located in relatively flat terrain with small hills and heterogeneous surface cover about 1 km southeast of the foot of Schetelig mountain (694 m MSL) and about 1.5 km from the foot of the Zeppelin mountain (556 m MSL). For more than 15 years, it has provided climate, ground temperature, and soil moisture data that have been the basis for studies, for example, on the ground thermal regime (Roth and Boike 2001; Boike et al. 2007; Westermann et al. 2010; Weismüller et al. 2011), the energy balance of the snow cover (Winther et al. 1999; Boike et al. 2003; Westermann et al. 2011), and the performance of satellite products (Westermann et al. 2010, 2011). Since 2007, the turbulent fluxes have been measured by the eddy covariance technique, as well as the fluxes of carbon dioxide (Lüers et al. 2014). A set of independent measurements of all components of the SEB for 1 year (Westermann et al. 2009) is used for model validation in this study. Measurements of incoming shortwave and longwave radiation (SWin and LWin) are taken from the Baseline Surface Radiation Network (BSRN) station in Ny-Ålesund (Ohmura et al. 1998). The outgoing longwave radiation (LWout) is measured at the Bayelva site, and the outgoing shortwave radiation (SWout) is taken from the BSRN measurements. The surface albedo is calculated using the BSRN shortwave radiation measurements (SWout/SWin), except during the snowmelt period when a constant surface albedo of 0.65 was used. Flux Qg is inferred from profile measurements of ground and snow temperatures. All of the components of the SEB from Bayelva are available at 1-h resolution.
Daily radiosoundings launched at the Alfred Wegener Institute for Polar and Marine Research/Institute Polaire Paul Emile Victor (AWIPEV) research base in Ny-Ålesund (Maturilli 2009) provide vertical temperature and humidity profiles for validation. Other measurements include the cloud-base height (CBH) retrieved by a laser ceilometer with 1-min time resolution (Maturilli 2011). Here, we use the median CBH value in the 15-min interval centered at each whole hour to compare with the instantaneous hourly output from WRF, using only the instances in which at least 80% of the measurements in this interval agree on the presence of clouds in a given part of the atmosphere. Measurements of precipitation, wind speed, 2-m air temperature T2, surface pressure, and relative humidity (RH) from the Norwegian Meteorological Institute (http://www.eklima.no) are used for further validation at this site.
The other set of validation data for SEB is from an automatic weather station (AWS) on the Austfonna ice cap (Schuler et al. 2014). The AWS is located at the western part of Austfonna (Etonbreen) at 370 m MSL. Here four components of radiation are measured, in addition to surface meteorological variables. Because of its remote location, the station is only maintained once per year during a field campaign in April/May, and therefore sporadic data gaps and episodes with lower data quality occurred. The AWS nevertheless gives valuable information about the SEB in a region where other observations are sparse.
3) Secondary sites
We also use data from a network of manual and automatic weather stations used by the Norwegian Meteorological Institute (http://www.eklima.no) to validate the general performance of the model. Here we use measurements of RH, T2, wind speed, and surface pressure, mainly to validate the model on the regional scale.
For our comparison, we first consider how well the model reproduces surface meteorological conditions over the entire archipelago (section 3a), showing results from both of the two inner domains (3 and 1 km) and using a constant lapse rate of −0.0065°C m−1 to correct for elevation differences between model and observations. For the detailed SEB comparisons at Bayelva (section 3b) and Austfonna (section 3c), only data from the finest available resolution are shown (1 and 3 km, respectively) with no height correction for temperature given that the elevation differences to observations are small (10 and 28 m, respectively). All results are from the nearest grid point with the right land cover in the model.
a. Surface meteorological conditions
The monthly T2 correlation coefficient exceeds 0.8 for all stations during most of the year (Fig. 2) but is lower during the three summer months [June–August (JJA) of 2008] and the first simulation month (15–31 March 2008). The temperature bias has the opposite form, with simulated mean values mostly within about 1°C of observations during the summer but otherwise approximately ranging from −3.5° to +3.5°C. As we will see in the results for Austfonna (section 3c and Fig. 9, described below), this follows the variation in temperature itself, with small variations during summer and larger ones during the rest of the year.
For the specific humidity q (calculated from observed RH) the correlation closely resembles the one found for T2, although with somewhat higher values (mostly >0.7). Also similar for T2 and q is an apparent gradient in the biases from mostly positive in the northeast to negative in the southwest.
The wind speed correlation is mostly between 0.4 and 0.8, with no clear seasonal or regional pattern, with monthly biases being almost exclusively positive, ranging from −0.6 m s−1 at Edgøya in October 2008 to 3.3 m s−1 at Hornsund in April 2008. Not shown in the figure is the surface pressure, for which the correlation is 0.99 for all stations and all months.
b. SEB at Bayelva
For the comparison of the SEB fluxes at the Bayelva site (Table 1), we divide the year into six seasons by following the method of Westermann et al. (2009). The summer season (1 July–31 August) is when the ground is (mostly) snow free and SW radiation dominates the SEB. The autumn season (1–30 September) is the transition season from SW-dominated SEB to the dark season and is often the time during which a snow layer is formed. The dark winter season (1 October–15 March) is dominated by LW radiation. During the light winter season (15 March–15 April), the SW radiation increases again, but its influence on the SEB is small because of high surface albedo and large zenith angles, and this is also the season with the lowest soil and skin temperatures Ts. Following the light winter period is the premelt season (16 April–31 May) during which the SW radiation becomes more important and the temperature increases. Last is the snowmelt season (1–30 June), during which the bare tundra starts appearing and the largest values of SWin occur. Although the different seasons have different characteristics for SEB, the transition between them is often gradual. More details on the different seasons can be found in Westermann et al. (2009).
In the following we use the convention that positive values denote fluxes away from the surface (to the atmosphere or to the soil) and negative values denote fluxes to the surface. Net SW and LW radiation (ΔSW and ΔLW) are given as the sum of the positive (outward) component and the negative (downward) component.
1) Summer (1 July–31 August 2008)
During the summer season, the observed temperatures are captured well by the model, with (negative) biases of less than 0.5°C for both T2 and Ts and correlations of 0.88 and 0.91 respectively (Table 1). Considerable differences are, however, found for the individual fluxes. ΔSW has a bias of −19 W m−2 but is partly compensated by a ΔLW bias of 8.3 W m−2. These biases are both mainly the result of distinct days that lack clouds in the model. Figure 3 shows the radiosonde profiles for the five days with the highest LWin biases (which mostly coincide with the days with maximum SWin biases), accounting for more than 40% and 55% of the LWin and SWin biases, respectively. On two days (5 July and 2 August), the model simulates a temperature inversion at approximately the right height, although weaker and with a too-warm PBL and no saturation. On 8 July and 19 August, the height of the temperature inversion is underestimated and the saturation seen in the upper part of the PBL is not captured by the model. On the final day (20 August), the model underestimates the humidity and overestimates the temperature in much of the lower atmosphere and misses a small inversion at ~2300 m (not shown).
The results also show a very high Bowen ratio (Qh/Qe, hereinafter BR) of ~6, as compared with ~1 in the observations. This result arises from a combination of a large overestimation of Qh (66 W m−2, as compared with 21 W m−2 in the observations) and an underestimation of Qe (12 W m−2, as compared with 23 W m−2 in the observations). Still, Qg is only overestimated by 1.3 W m−2.
Model adjustments: Soil hydrology for frozen ground
To understand the strong overestimation of BR, two sensitivity simulations were performed for the summer season with adjustments in the soil-hydrology calculations in the Noah LSM scheme. These were both started 8 days before the summer season (i.e., 23 June) to include the effect of water from snowmelt.
The first adjustment was made to the calculation of surface runoff, which originally was reduced to account for frozen water in the soil column without taking into account which soil layers were frozen (Koren et al. 1999). For the Bayelva site where ice remained in the lowest soil layer throughout the summer, this resulted in about 80% of the precipitation during this season being removed as surface runoff, even though the upper layers of the model were nonfrozen and relatively dry. In the sensitivity test (referred to as “SURFADJ”), this effect was removed so that the model no longer reduced infiltration when the ground was frozen.
Second, an adjustment to the treatment of underground runoff was made so that no water could be removed through the bottom layer when it contained more than 50 mm of frozen water. In addition, the model was modified to allow for more water infiltration at the surface (minor adjustment). The sensitivity simulation called “FULLADJ” includes all of these modifications (Figs. 4 and 5).
These modifications drastically increased the soil water content and decreased the simulated BR (Fig. 4). Removing the frozen-ground infiltration reduction (SURFADJ) made a large difference during the rain events around 10–13 July and 10–15 August. For the underground-runoff adjustment, the difference is most clear during periods without precipitation from 1 to 9 July and from 14 to 18 July, when the two modified runs lose water at about the same rates but with considerably lower BR for the simulation without underground runoff (FULLADJ). From about 25 July, the soil becomes relatively dry again in all three simulations and the BR increases to above 5, even with the observed value staying around 2. From the second large rain event (10–15 August), the BR in both adjusted simulations drops to values that are close to the observed values of ~1, although with a consistent (small) positive bias in the SURFADJ simulation.
The seasonally averaged energy fluxes also show decreasing BR with each modification, with the FULLADJ simulation resulting in a BR of ~2 (Fig. 5). This result is closer to the observed value of ~1. With the precipitation being overestimated in this season (Table 1), however, one would expect the simulated BR to be too low rather than too high. Also clear from Fig. 5 is the very small sensitivity of the other SEB fluxes to these modifications, with net radiation only increasing slightly and Qg staying approximately the same.
2) Autumn (1–30 September 2008)
Very high correlations are found for T2 and Ts during the autumn season (0.93 and 0.98, respectively), although both show negative biases of −0.8°C (Table 1). The radiative fluxes also have high correlations (above 0.8) for all components except the SWout.
The low correlation and relatively high bias in SWout indicate a wrong timing of the forming of a snow layer. In the simulations, the snow-covered fraction increases from 0.0 to ~0.4 on 23–24 September, with a corresponding increase in surface albedo from 0.15 to 0.43. This increase in albedo is not found in the observations.
The simulated soil moisture in the first half of this season is too dry (not shown), with similar underestimation as seen in the last part of the summer (“ORIGINAL” in Fig. 4b). This results in the simulated Qe being only ~40% of the observed value in this season. It also affects Qh, which is too high in the beginning of the month but is compensated by a too-strong negative flux (to the ground) in the last part of the month, giving a net negative bias.
3) Dark winter (1 October 2008–15 March 2009)
For the period of dark winter we again see good agreement for temperature, with T2 having a correlation of 0.95 and a bias of 0.6°C (Table 1). For Ts, the correlation is 0.91 and the bias is −1.2°C, which results in a negative LWout bias of −6.1 W m−2. The underestimation of Ts is in line with a positive bias (underestimation) in LWin of 14 W m−2.
To further investigate the bias in LWin, the hourly data are studied in relation to the existence of low clouds (<3 km) in the model and in ceilometer data, defining clouds in the model as any layer with cloud water content (CWC) of more than 10 mg m−3 (Fig. 6). This reveals a general underestimation of cloudy hours. During about 22% of the dark winter season, clouds can be found only in the ceilometer data (black dots), whereas clouds are found only in the model approximately 4% of the time. During most of the time when the model incorrectly simulates no low clouds, saturation is reached also in some model layer below 3 km, but with too little CWC to have a significant radiative effect. For the correctly simulated low clouds, the hourly data show a substantial spread, with a mean LWin bias of 4.1 W m−2. A certain spread is also found for the correctly simulated cloud-free conditions (blue dots), with a mean LWin bias of 2.2 W m−2.
Larger simulated than observed magnitude of Qh and Qe is also found in this season (Table 1). Even though they are of different sign, and should therefore have opposite dependence on the Ts, the simulated values of −37 and 6.8 W m−2 are both larger by more than a factor 2 than the observed ones, with correlations of only 0.36 and 0.48, respectively.
The simulated Qg (−4.0 W m−2) is close to the average measured snow heat flux (−5.0 W m−2) in this season, despite the simple treatment of snow in Noah and the observations giving the snow heat flux in only the upper part of the snowpack.
4) Light winter (15 March–15 April 2008)
The largest temperature biases are found in the period of light winter at 3.0° and 5.0°C for T2 and Ts, respectively (Table 1). The bias in Ts, combined with a negative bias in LWin, gives rise to a ΔLW that is more than 2 times the observed value. We note here that the observational data from Bayelva during this season have considerable gaps, however, and deviate substantially from the observations from the maintained BSRN station and should therefore be considered to be very uncertain.
Flux Qh is the main component compensating for the radiation loss and is again overestimated relative to the observations. The magnitude of Qe is also overestimated (by a factor of 4), although it is still a small component in the SEB. Also in line with a too-high Ts is the underestimation of Qg that is seen in this period.
5) Premelt (16 March–31 May 2008)
During the premelt period, we find the largest biases in the individual radiation components, in both absolute and relative values (Table 1). The SWin is overestimated by 56 W m−2 (30%) relative to the observations. The effect of this high SWin bias is substantially suppressed by the high surface albedo, even though the simulated surface albedo is slightly lower than the observed (0.72 vs 0.78). Together this gives a ΔSW bias of −26 W m−2, which is almost entirely balanced by a positive ΔLW bias of 22 W m−2, both of which are consistent with too few clouds or clouds that are too optically thin.
For the other SEB fluxes, only Qe shows a large bias in this season (14 W m−2), whereas the biases in Qh and Qg are only 1.1 and −2.6 W m−2, respectively. Still both turbulent fluxes show considerable mean absolute error (MAE). To understand these differences better, we look more closely at the first 15 days of this season (Fig. 7). These days can be divided into a relatively cloudy period (16–23 April) with measured 10-m wind speeds mostly above 6 m s−1 followed by a consistently calm period (wind speeds < 5 m s−1) with mostly clear skies (24–30 April). During the first part the simulated Ts is in fair agreement with observations, and both the simulation and the observations show low values of Qg, but the model overestimates the magnitude of both Qh and Qe. In the second part we see, however, that Qg becomes an important part of the observed SEB, whereas it is still close to zero in the model. At the same time, the diurnal cycle of Ts is considerably amplified. In a similar way, Qh and Qe are both overestimated during daytime, but because Qh is compensated with negative values during night, a large bias is only seen for Qe.
6) Snowmelt (1–30 June 2008)
In the snowmelt season there are also radiation biases that are consistent with too few clouds or optically thin clouds, with average biases of −56 and 13 W m−2 for SWin and LWin, respectively (Table 1). The net radiation bias is only −6 W m−2, however. As in previous seasons, the magnitudes of both turbulent fluxes are overestimated in this season, but to some extent these biases cancel out because they have different signs.
A major factor controlling the SEB is the timing of the snowmelt, mainly because of the large associated change in surface albedo. This factor is especially important in this region, because it coincides with maximum SWin around the summer solstice. Figure 8 shows that the timing of the snowmelt is captured remarkably well, as seen both from when the snow height becomes zero and from when Ts becomes (significantly) positive. Comparison between simulated and measured amounts of snow melted in this season suggests that this is partly a result of compensating errors, however. From the average snow depth and seven snow-density measurements in the study area, Westermann et al. (2009) estimated the average snow water equivalent (SWE) shortly before the snowmelt season (25–28 May) to be ~210 mm, as compared with 310 mm in the model. If all of this SWE is melted during the snowmelt period, it should give an average energy consumption by snowmelt of ~39 W m−2 (which is close to the simulated residual), as compared with 27 W m−2 with the measured amount of snow. This difference in SWE seems to originate from the initialization with the coarse-resolution spinup. Comparison with the nearby grid points with snow depths similar to the observed one shows a transition to mostly snow-free ground 5–6 days earlier than what is seen in Fig. 8, with a considerable effect on the SEB.
7) The annual SEB
Annually averaged, both the temperatures and the individual SEB fluxes are captured well by the model (Table 2), despite the considerable biases in several seasons. The simulated annual average T2 (−4.9°C) and Ts (−7.4°C) are both within 0.5°C from the observed values and have MAE of 2.0° and 2.9°C, respectively.
The model correctly simulates ΔSW (−46 W m−2) as the dominant source of energy to the surface and ΔLW (43 W m−2) as the dominant sink. These are both ~20%–25% larger than observed on average, but together they give a net radiation bias of only 2 W m−2.
As in the observations, the simulated turbulent fluxes are also of opposite sign and add up to nearly zero, although turbulent fluxes are larger in magnitude in the simulation. The relatively small annual biases in these fluxes are, however, a result of compensating errors giving large MAE of both. The small annual bias in Qh is a result of it changing sign between summer and winter, whereas for Qe the underestimation during summer and autumn that is due to the too-dry soil compensates for the overestimation throughout the rest of the year. The annually averaged Qg is small in both the simulation (1.1 W m−2) and the observations (0.8 W m−2), which is what one would expect for a permanent permafrost site like Bayelva.
c. SEB at Austfonna
At Austfonna the T2 correlation is high (~0.8 during most months) and the biases are mostly less than ±1°C (Fig. 2). During the dark season (last 6 months of the simulation), however, the bias increases to between −2° and −4°C, even though the correlation stays high. Figures 9a,b reveal that this period is characterized by large temperature variations, with both T2 and Ts changing by more than 20°C within days. These rapid changes are captured well by the model, although the temperatures are generally underestimated in this season.
An overestimation of the annual temperature cycle is seen, with especially Ts being too high during the summer and too low during the winter. This corresponds with the bias in net radiation (Fig. 9c), with too-strong net downward (negative) flux during summer and too-strong net upward (positive) flux during winter. As was the case for Bayelva, this is consistent with too few clouds or with clouds that are too optically thin. In addition, the model generally underestimates the surface albedo at this site (Fig. 10), which contributes to the radiation bias during spring and summer, even though the variation in surface albedo is to some extent captured.
The simulated temperatures (both Ts and T2) were found to be in good agreement with observations, especially at the two main locations (Figs. 2 and 9 and Table 1). At Bayelva, the temperature biases in the individual seasons were mostly ±1°C or less and correlations were never below 0.88. With updated high-quality SST and sea ice fields, the model therefore seems to downscale the temperature from the reanalysis in this region well, even without nudging or reinitialization.
For the individual SEB fluxes, substantial biases and MAE were found in some seasons. In addition, the MAE of Ts was considerable in parts of the year, especially the coldest seasons. In the following, we will discuss these differences in relation to how the model simulates clouds, surface albedo, soil moisture, and turbulent fluxes before finally discussing the issue of spatial variability within the study area.
Several studies have shown that reanalyses have problems reproducing observed cloud conditions in the Arctic (Walsh et al. 2009; Chernokulsky and Mokhov 2012; Zib et al. 2012), with even the annual cycle of the cloud fraction differing among individual reanalyses. Ground-based observations consistently show a maximum cloud amount during the “Arctic stratus” season in summer and early autumn (Warren et al. 1988; Intrieri et al. 2002). Arctic stratus clouds are a challenge for atmospheric models because their persistence relies on a complex interaction between several processes: boundary layer turbulence, cloud-top radiative cooling, cloud microphysics, and humidity advection (Sedlar and Tjernström 2009). In addition, the Arctic low clouds are very often of mixed phase, with a delicate balance between supercooled water and ice crystals (Morrison et al. 2012). The models often fail to account properly for this complex interplay of processes, and, as a result, key cloud characteristics such as liquid water path, ice fraction, and cloud droplet size often have significant biases (Morrison et al. 2009). In their simulation of the Svalbard region with the WRF Model, Claremar et al. (2012) found considerable radiation biases that could be related to clouds and also showed that this result was highly dependent on model resolution. Our results therefore seem to be as good as could be expected with current regional or local-area models. Still, because clouds have a dominating effect on the SEB in the Arctic (Langer et al. 2011a,b), it is clear that simulating these clouds remains a major challenge for SEB calculations in this region, with our analysis showing that there are different problems in the different seasons.
During summer, the radiation biases in Ny-Ålesund were found primarily to be a result of distinct days with incorrectly simulated cloud-free conditions in the model. Common features of the five days with the largest LWin biases were differences in the temperature and/or humidity profiles in the PBL between model and observations (Fig. 3). The model also underestimated the amount of time with low clouds present during the dark winter period (Fig. 6), although during this season tiny amounts of CWC were still present most of the time when clouds were observed. In addition to fewer cloudy hours in the model, a mean underestimation of LWin was found during the correctly simulated occurrences of low clouds and to a lesser extent also during the correctly simulated cloud-free conditions. A more detailed analysis of this bias would require looking closer at cloud lifetime, total water path, ice fraction, and droplet size, which is beyond the scope of this study. The bias during cloud-free conditions cannot be ignored, although its contribution to the mean SEB in the different seasons appears to be less important than the simulated cloud properties.
b. Surface albedo
During most seasons the average surface albedo at Bayelva is captured well, and at Austfonna its evolution is qualitatively in agreement with observations. The highest values of snow albedo during the spring are substantially underestimated at both sites, however, showing that 1) the maximum fresh-snow albedo of 0.85 is too low in this environment, 2) the land-cover-dependent snow albedos are too low for the respective land types, and/or 3) the weighting between these two (here equal) should be more toward the former. Long records of measurements of surface albedo from both Austfonna (Schuler et al. 2014) and Ny-Ålesund (Winther et al. 2002) show that values above 0.8 are often observed at both sites during spring and that values above 0.85 are not uncommon at the glacier site. In addition, the lowest simulated values of 0.7 for ice surfaces (seen in Fig. 10) are usually not low enough at Austfonna, where values down to about 0.4 are found when all of the snow is melted (Schuler et al. 2014). In part, this high minimum value is because the Noah scheme uses a fixed lower value for SWE on glaciers (Collier et al. 2013), and hence bare ice is never really exposed. With an adjustment of the maximum value and removal of the limit to the minimum SWE, we consider this albedo scheme to be sufficient to reproduce seasonally averaged values, although a more sophisticated scheme would be needed to capture changes on shorter time scales or responses to, for instance, changing wind and temperature conditions.
c. Soil moisture
From the sensitivity studies of soil hydrology during summer, it is clear that the surface and underground runoff dominated the soil hydrology in the original model, in clear contrast to the observations in which the summer evaporation exceeded the precipitation. The first model adjustment—removal of the reduction in surface infiltration that is due to frozen ground—has a clear physical basis because frozen water below 1-m depth cannot prevent the water from penetrating into the top soil layers. Our solution of entirely turning off this effect could lead to an overestimation of infiltration when the first soil layer is actually frozen. In the future, the existing formulation should therefore be adjusted to account for which part of the soil column is frozen.
For the underground runoff, there are physical processes that could result in this effect, even when there is permafrost present. These processes include horizontal groundwater flow along the top of the permafrost in sloping terrain or penetration into the permafrost if it is discontinuous or inhomogeneous. Our approach of removing underground runoff when frozen water is present can therefore not always be justified. At Bayelva, however, the loss of soil water through underground runoff should not be a major contribution on the scales used in our simulations, although the effect of sloping terrain can have an effect near the mountains and at subgrid scales.
d. Turbulent fluxes
The lowest correlation and largest deviations from observations are found for the turbulent fluxes, both annually averaged and within most seasons. Some of these deviations must be considered to be a result of deficiencies in other parts of the model than the PBL scheme, in particular the radiation bias and the soil moisture deficit discussed above. In addition, the observed residual is considerable in many seasons and follows the sign of the dominant turbulent fluxes during much of the year (e.g., summer and dark winter). Foken (2008) described how a residual is typically found when measuring SEB in heterogeneous landscapes and attributed it to secondary eddies contributing to Qh and Qe that cannot be measured with the standard eddy covariance technique. We must therefore anticipate that the values of these fluxes in the observations in several seasons are too low. The differences between simulated and observed turbulent fluxes are hence related not only to how these fluxes are simulated in the model, especially during summer. Still, the large MAE, relatively low correlation, and persistent overestimation of both Qh and Qe reveal that there are real problems with simulating these fluxes, especially in the stable boundary layer (SBL) that is typical for most other seasons in this area.
The simulation of the SBL is a problematic issue in both weather and climate modeling (e.g., Holtslag et al. 2013). Considerable effort has been put forth to improve this, for instance, through the different GEWEX Atmospheric Boundary Layer Study (GABELS) experiments (Holtslag 2006; Holtslag et al. 2012). Operational as well as research models differ in both how the equations for turbulence are closed and the stability functions used to account for atmospheric stability (e.g., Cuxart et al. 2006), and many models are known to overestimate turbulent fluxes in the SBL relative to measurements.
For regional simulation of Svalbard, the sensitivity to selected PBL parameterizations that are available within the WRF Model has been assessed by both Claremar et al. (2012) and Kilpeläinen et al. (2012), with the former study finding mean T2 differences of up to 1°C during winter. In a more detailed study of the SBL with a single-column model, Sterk et al. (2013) found different processes to be important in different wind regimes. During high wind speeds, the simulated SBL was found to be most sensitive to the mixing processes (calculated by the PBL and surface-layer schemes), whereas coupling to the surface (through Qg) and radiation were found to be more important in calm conditions.
Our results also indicate different sources of error in the SBL during different wind regimes. During the calm, clear-sky period at the end of April (Fig. 7), the overestimation of diurnal cycle of both Ts and turbulent fluxes must at least partly be attributed to underestimation of Qg variations owing to the simple snow treatment. With too much energy at the surface during daytime, the model simulates too-high Ts and too-large turbulent fluxes. During night this is to some extent compensated for by underestimation of Ts and too-large negative values of Qh, giving relatively small biases but considerable MAE for these variables. On the other hand, Qe is not compensated for by negative values during night and therefore shows both a large bias and a large MAE in this season. Underestimation of Qg in the presence of thick snow in Noah was also noted by Niu et al. (2011) and was part of their motivation for developing the new Noah Land Surface Model with Multiparameterization Options (Noah-MP) with several separate snow layers.
The overestimation of turbulent fluxes is a persistent feature, however, and not only when Qg is underestimated, and therefore at least two more factors need to be considered here. First, these fluxes are sensitive to the values of roughness lengths (see section 2a). Whereas summer measurements of Z0m in this region are mostly smaller than our value (Lloyd et al. 2001; Sjöblom 2014), there is support in the literature for using a smaller value of Z0t relative to Z0m in this kind of environment (Chen and Zhang 2009). The direct effect of decreasing (increasing) either is smaller (larger) values of Qh and Qe, but the full effect can only be assessed through a full simulation because it influences the whole structure of the PBL. Second, we emphasize that our simulations are performed with relatively coarse vertical resolution, which could be a problem when simulating SBLs (e.g., Steeneveld et al. 2006). In particular, the very shallow inversion layer close to the ground that has been measured in this region during spring cannot be expected to be adequately represented in our simulations and might even influence the measurement of Qh and Qe (Lüers and Bareiss 2010).
e. Spatial variability
Even with horizontal grid spacing in the model of 1 km, we also need to consider the spatial variability on the subgrid scale, and to what extent the SEB measurements (mainly point measurements) are representative for such an area. By mapping the snow depth at Bayelva with ground-penetrating radar and manual-probe measurements, Gisnås et al. (2014) found the snow depth to vary from almost bare blown to 2 m within an area of about 1 km2, with the corresponding mean annual ground surface temperature varying between −5° and 0°C. In a similar way, Westermann et al. (2010) found that the thaw depth varied between 1.6 and 2.0 m along a 175-m transect in the same area at the end of the summer season in 2008. Together these results illustrate that the SEB is a local feature, which must be kept in mind when comparing measurements and model output even at a relatively fine resolution.
In this study, we have compared WRF simulations of all components of the surface energy balance at Svalbard with in situ measurements. The model was found to reproduce monthly-to-seasonal mean temperatures throughout this region very well, and the characteristics of the SEB in the different seasons were also broadly in agreement with observations. Still, several issues were found in the simulation of the individual energy fluxes:
During much of the year, large biases were found in radiation components, mainly because of too few clouds or clouds that are too optically thin in the model. During summer, these biases could mainly be related to distinct days on which the model erroneously simulated cloud-free conditions and observations showed clouds in the upper part of the PBL, whereas the model seemed to simulate winter clouds as too optically thin.
During summer, the model initially simulated a far too-high Bowen ratio because it simulated too little soil water. Two physically based modifications to the soil hydrology for frozen ground improved the simulations greatly in this respect.
Our results suggest that underestimation of the magnitude of ground heat flux Qg on shorter time scales (e.g., the diurnal variation) resulting from the crude snow treatment leads to overestimation of diurnal cycles of skin temperature and can give biases in sensible and latent heat fluxes, even when Qg itself does not show a significant bias over the length of a season.
Too-low simulated snow albedos were found to enhance the net shortwave radiation bias, especially at the Austfonna ice cap during spring, showing that higher maximum snow albedos should be used in this environment.
We therefore conclude that, although it reproduces the mean temperatures well in this region, there are limitations to using the model directly for permafrost and glacier modeling because of the considerable biases in individual energy fluxes. Some weaknesses could be identified and corrected with the kind of model evaluation done here, but adequately simulating Arctic clouds and adequately simulating turbulent fluxes in the stable boundary layer remain the largest unresolved issues in the quest to correctly reproduce the SEB in this region.
This research was done as a part of the CryoMet project funded by the Norwegian Research Council (NFR 214465). In addition, a field trip to Longyearbyen and Ny-Ålesund was financed by Svalbard Science Forum (through the CRYO-FImBack workshop 2013; NFR 226508) and Department of Geosciences, UiO, respectively. The authors are grateful to three anonymous reviewers for valuable comments and suggestions.
List of Symbols
SWin Incoming shortwave radiation
SWout Outgoing shortwave radiation
ΔSW Net shortwave radiation
LWin Incoming longwave radiation
LWout Outgoing longwave radiation
ΔLW Net longwave radiation
Qh Sensible heat flux
Qe Latent heat flux
Qg Ground heat flux
T2 2-m air temperature
Ts Skin temperature
q Specific humidity
Res Residual of surface energy balance
The choice of changing radiation schemes between the inner two domains might not be optimal, but sensitivity simulations indicate that it is not important for the results.