2.1 Downscaling

Precipitation is often heavily biased in coarsely resolved reanalyses, especially in environments with pronounced topography, where it typically is too low and lacks spatial detail (Schuler et al., 2008). This is associated with the smoothed representation of the actual topography in the large-scale model used for the reanalysis (Fig. 1), leading to an underestimation of orographic precipitation. Figure 1 shows that ERA greatly generalizes the high-resolution topography, representing Svalbard as a wide and flat bump that exceeds sea-level far off the actual coastlines, while surface elevation in the interior does not exceed 400 m a.s.l. In contrast, the highest elevation in our gridded topography map is 1600 m a.s.l. The roughness of the actual topography that gave rise to the name of the main island “Spitsbergen” (pointed mountains) is not represented by the smoothed topography used for the ERA reanalyses. We assume that this is the main reason for the poor performance of reanalyzed precipitation. To account for orographic enhancement, we use a linear theory (LT) of orographic precipitation (Smith and Barstad, 2004).

The other required climate variables are downscaled to the 1 km grid largely following the TopoSCALE methodology (Fiddes and Gruber, 2014), which also builds on the assumption that weaknesses in the representation of topography at the coarse scale are mainly responsible for the misfit between coarse-scale and point observations. TopoSCALE exploits the relatively high vertical resolution of the reanalysis data to downscale variables to the elevation of the actual topography, based on the properties of the vertical structure in the reanalysis. The downscaled fields preserve the horizontal gradients present in ERA but include additional features caused by the real topography that were not present in the ERA products. In doing so, we add spatial detail to the reanalysis fields that is consistent with the temporal evolution of atmospheric conditions of the reanalysis. This approach is assumed to outperform simpler bias corrections, since transient properties of the atmosphere are accounted for. For example, transient lapse rates including inversions in the reanalysis data will be preserved in the downscaled product.

In our application, we modified the TopoSCALE methodology regarding downscaling of direct shortwave radiation and air temperature, as described in Sect. 2.3 and 2.4.

2.2 Precipitation

The LT model describes an air parcel as it moves across a prescribed surface topography. The air parcel is characterized by its temperature, stability, wind direction and speed. Terrain-induced uplift of the air parcel results in condensation and eventually precipitation of moisture downstream of the uplift. This model has been successfully evaluated using precipitation gauges (Barstad and Smith, 2005) and snow measurements (Schuler et al., 2008; Østby et al., 2017) and applied for downscaling precipitation (e.g., Crochet et al., 2007; Jarosch et al., 2012; Roth et al., 2018). The linear theory utilizes a Boussinesq description of mountain wave to derive a transfer function that, for given wind conditions, relates the orographically enhanced precipitation to terrain topography.

By using spectral decomposition and algebraic manipulation, Smith and Barstad (2004) derived the following transfer function:

relating $\widehat{\mathbf{P}}(k,l)$, the Fourier transform of the precipitation enhancement, to the Fourier transform of terrain elevation $\widehat{\mathbf{h}}(k,l)$, with k and l being the horizontal wave numbers. This relation depends on the uplift sensitivity factor C_w, thickness of the moist layer H_w, the intrinsic frequency $\mathit{\sigma }(k,l)=Uk+Vl$ (U and V being the east and north components of the wind vector), and the conversion and fallout timescales τ_c and τ_f, respectively. In Eq. (1), the vertical wave number m controls the depth and tilt of the forced air uplift and is a function of the moist Brunt–Väisälä frequency, N_m, a quantity describing atmospheric stability.

Precipitation rates are obtained by retransforming $\widehat{\mathbf{P}}$ and adding it to the background precipitation P_∞, which accounts for large-scale frontal and convective precipitation separate from orographic precipitation P_oro. P_∞ has been corrected for the orographic effect already present in the ERA reanalyses P_oro(h_ERA) by estimating this effect applying Eq. (1) to the large-scale topography h_ERA and removing the result from the ERA precipitation ${\mathbf{P}}_{\mathrm{\infty }}={\mathbf{P}}_{\mathrm{ERA}}-{\mathbf{P}}_{\mathrm{oro}}\left({\mathbf{h}}_{\mathrm{ERA}}\right)$ (Schuler et al., 2008). Total precipitation, P_total, is then

Since the theory assumes saturated conditions, we account for reduced orographic enhancement at lower humidity by adopting a correction factor f proposed by Sinclair (1994):

which suppresses orographic enhancement when RH <0.8.

Instead of treating N_m, τ_f and τ_c as adjustable, constant parameters, we exploit the evolution of the moisture-bearing layer of the atmosphere described in the reanalyses to derive transient values. In doing so, we remove calibration parameters from our method and enable weather-dependent variation in N_m, τ_f and τ_c. Values of N_m are calculated as follows:

where g is gravitational acceleration and $\stackrel{\mathrm{‾}}{T}$ is vertically averaged air temperature weighted by the moisture content at several pressure levels (Jarosch et al., 2012). Environmental lapse rates Γ_e are derived from air temperature at 700 and 850 hPa and corresponding geopotential heights. Γ_m is the moist adiabatic lapse rate, calculated according to Stone and Carlson (1979) using vertically averaged values of atmospheric properties from the reanalyses weighted by moisture content. This follows the convention that a positive lapse rate represents cooling with increasing elevation. Barstad and Smith (2005) report that typical values of N_m range between 0 s⁻¹, representing an atmosphere with no stratification, and 0.01 s⁻¹, representing a stably stratified atmosphere. To avoid conditions inconsistent with the assumptions of the theory, we limit N_m to this range. The quantities H_w, C_w and m are derived from N_m (Smith and Barstad, 2004).

Advection timescales τ_c and τ_f are assumed equal, and $\mathit{\tau }={\mathit{\tau }}_{\mathrm{c}}={\mathit{\tau }}_{\mathrm{f}}=H_w/v$ is derived from the thickness of the moist layer H_w and accounting for a typical hydrometeor fall speed v, which is taken as constant but allowed to take different values for solid and liquid hydrometeors. This phase transition is determined by a threshold temperature of 273 K, hence v(T≤273 K) =1 m s⁻¹ for solid and v(T>=273 K) =2 m s⁻¹ for liquid precipitation. In Eq. (1), terrain elevation h is the only gridded variable, the other variables represent averages over the volume of the described air parcel. To characterize this air parcel, we first vertically average the values defined at the nodes of the horizontal domain over the 700 and 850 hPa pressure levels weighted by the moisture content of the individual layers. These vertically averaged values are then horizontally averaged over an area defined by a 200 km buffer around the 200 m contour of the reanalysis topography, the latter roughly outlining the extent of the archipelago (Fig. 1).

This setup is then applied to each 6 h time step, and the resulting time slices are progressively added to the record.

Figure 2Example plot of long-term mean annual precipitation (m) over 1979–2017, downscaled from ERA-Interim.

2.3 Air temperature

The downscaling for near-surface air temperature at the 2 m level (T_2 m) closely follows the TopoSCALE procedure of Fiddes and Gruber (2014), thereby we assume that the vertical structure of the free atmosphere determines the distribution of T_2 m with terrain elevation. For each 6 h time step, T_2 m is derived from a three-dimensional interpolation of the vertical air temperature structure of the large-scale reanalysis to the location of the grid nodes representing the high-resolution terrain elevation. We notice that for a melting snow or ice surface, skin temperature is bounded to 273 K, influencing the near surface air temperature and resulting in reduced along-surface air temperature lapse rates (e.g., Marshall et al., 2007). To account for this effect, we apply a simple horizontal interpolation of the two-dimensional T_2 m from the ERA reanalysis field, where T_2 m>273 K, instead of interpolating a three-dimensional data volume to the surface elevation of the high-resolution topography. This strategy is motivated by the discovery of unrealistic warm air temperatures at higher elevations in a test application. The occurrence of this unphysical temperature inversion was restricted to the melting period, caused by extrapolation of surface inversions close to a snow and/or ice surface, the temperature of which is capped at 273 K. To avoid this effect, we assume that where T_2 m>273 K, the T_2 m of the reanalysis is consistent with a melting surface and hence more realistic than a free-atmosphere interpolation.

Figure 3Example plot of mean air temperature (2 m) in ^∘C, over the period 1979–2017, downscaled from ERA-Interim.

2.4 Radiation

Downscaling of shortwave and longwave downwelling radiative fluxes (SW and LW, respectively) was conducted by adopting the TopoSCALE methodology (Fiddes and Gruber, 2014) with a few adjustments. The shortwave radiation at the surface level of the reanalysis is projected to the high-resolution topography in a three-step procedure: first the surface SW flux is separated into direct and diffuse components; second, the direct component is corrected for the elevation difference between the reanalysis surface and the high-resolution topography, considering an effective atmospheric transmissivity that is derived from top-of-atmosphere and surface fluxes; third, a topographic correction is applied to account for effects of slope and aspect of the high-resolution topography, as well as shading by surrounding topography. To compute direct solar radiation we apply the relationship of Kumar et al. (1997) for atmospheric attenuation rather than the one given by Fiddes and Gruber (2014). Solar geometry variables such as solar zenith and azimuth and topographic shading due to local slope and aspect are calculated following Reda and Andreas (2004). Cast shadow and hemispherical obstructions caused by surrounding topography are calculated following Ratti (2001). The longwave surface flux is downscaled by correcting for the elevation difference between reanalysis and high-resolution grids using an atmospheric emissivity. This emissivity is estimated by accounting for a clear-sky component that depends on humidity and air temperature and a cloud component that is estimated from the difference between the clear-sky component and the reanalysis longwave flux. Further terrain effects are incorporated through multiplication with the sky-view factor.

Figure 4(a) Example plot of 1979–2017 incident shortwave radiation (W m⁻²), downscaled from ERA-Interim. (b) Comparison between measured (Maturilli et al., 2015) and downscaled daily values of SW at Ny-Ålesund. The color indicates the density of points from red (high) to green (low).

Figure 5(a) Example plot of 1979–2017 downwelling longwave radiation (W m⁻²), downscaled from ERA-Interim. (b) Comparison between measured (Maturilli et al., 2015) and downscaled daily values of LW at Ny-Ålesund. The color indicates the density of points from red (high) to green (low).

2.5 Relative humidity and wind speed

Similar to our assumption about T_2 m over a melting surface, we suggest that wind speed and RH in the boundary layer are more affected by surface rather than by free-atmosphere conditions, and we hence apply a simple two-dimensional interpolation of the near-surface values of the reanalysis, instead of interpolating a three-dimensional data volume to the surface elevation of the high-resolution topography.

3.1 Precipitation

At the operational weather stations (Fig. 1), downscaled precipitation (Fig. 2) is overestimated by 5 to 25 mm per month at the weather stations, with a slightly higher bias during winter (Table 1). This is also consistent with the findings of Vikhamar-Schuler et al. (2019), who evaluated the performance at six weather stations for several 30-year reference periods (1961–1990, 1971–2000 and 1988–2017). These biases are partly caused by precipitation measurements that are too low, the values of which are heavily affected by wind-induced undercatch, especially for solid precipitation. Førland and Hanssen-Bauer (2000) suggest that for solid precipitation, the actual precipitation at wind-exposed sites may be up to 80 % higher than the gauge record. A newly developed correction scheme for Norwegian mountain environments (Wolff et al., 2015) supports the finding of a large undercatch for solid precipitation; however, this correction has not yet been applied for Arctic conditions in Svalbard.

In addition to gauge measurements from low-elevation stations along the coast, we also used snow survey transects across Austfonna, a large ice cap in northeastern Svalbard (Fig. 1), to evaluate the Sval_Imp precipitation. We suggest that snow deposition on large glacier areas, measured at the end of the winter, represents seasonally integrated precipitation. While these measurements do not allow temporal resolution below one snow season (typically October–May), they provide useful information about spatial precipitation patterns in areas and elevations not covered by the operational meteorological stations. This approach builds on the implicit assumption of negligible sublimation, such that accumulated snow water equivalent represents the sum of precipitation over the winter season. Svalbard has high air humidity throughout the year, limiting the potential for sublimation. In an energy-balance study, Østby et al. (2017) estimated sublimation to about 0.016 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$, which is 1–2 orders of magnitude smaller than typical annual precipitation sums. There is generally good agreement concerning the spatial pattern across the Austfonna ice cap (Fig. 6), where snow accumulation reveals a distinctive southeast–northwest asymmetry (Taurisano et al., 2007; Schuler et al., 2007; Dunse et al., 2009) caused by orographic enhancement of precipitation coming from the southeast sector (i.e., the Barents Sea), although in individual years the downscaled values underestimate the measured accumulation, especially in 2007 (Fig. 6). Even though there is considerable scatter between observed and downscaled winter precipitation, there is positive correlation, indicating that the spatial pattern is matched and in most years there is no systematic bias, showing that the overall precipitation amount is adequately represented. On the other hand, the unscaled ERA-Interim winter precipitation shows almost no spatial variation and considerably underestimates observed values (Fig. 2).

Østby et al. (2017) found that the winter mass balance of Hansbreen, a glacier close to Hornsund, was not well reproduced, both in terms of spatial pattern and accumulation amount. Aas et al. (2016) similarly reported the lowest performance for Hansbreen, although they used a much more complex precipitation scheme than the one presented here. The generally low performance of several precipitation distribution schemes compared to the Hansbreen record has been interpreted to result from local conditions at Hansbreen, where the spatial distribution of snow is caused by wind redistribution rather than by the spatial precipitation pattern (Grabiec et al., 2006).

Figure 6(a) Map view of Sval_Imp precipitation accumulated over October 1998 to March 1999, overlaid with colored circles that indicate the measured snow water equivalent (SWE) by Sand et al. (2003). (b) Scatterplot comparing the 1999 measurements to winter precipitation according to Sval_Imp (crosses) and ERA-Interim (dots). (c) A similar scatterplot to (b) but only for the Austfonna dataset, which provides multi-temporal coverage (Taurisano et al., 2007; Dunse et al., 2009).

3.2 Air temperature

Downscaled air temperatures (Fig. 3) are compared to observations at the meteorological stations listed in Table 1 mostly for the period after 2004. Despite altitude differences of up to 100 m between measuring site and corresponding grid node in the model, no altitude correction is performed, due to unknown lapse rates. In general the agreement is good between downscaled ERA and observed air temperatures, with biases mostly below 1.5 K (Table 1). Despite a small bias for mean annual temperatures, there is a clear seasonal bias, with ERA temperatures too warm during winter and too cold during summer (Table 1). Although the biases are negative during summer, ERA is too warm over the glaciers during summer, when 2 m air temperatures are above freezing. These findings are consistent with those of Vikhamar-Schuler et al. (2019), who evaluated differences in seasonal mean values for different 30-year periods 1961–1990, 1971–2000 and 1988–2017.

At Svalbard Airport, the performances of downscaled ERA-40 and ERA-Interim are investigated for the entire model period. Over 1957–1979 only monthly measured air temperatures are available at Svalbard Airport, where downscaled ERA-40 has a monthly root-mean-square error (RMSE) of 2.3 ^∘C. For the 1979–2002 period the reanalysis products overlap with monthly RMSE of 1.8 and 1.5 ^∘C at Svalbard Airport for ERA-40 and ERA-Interim, respectively. We attribute the lower performance prior to 1979 to the lack of satellite observations to constrain sea surface temperatures and sea ice cover in the reanalysis. Since the Svalbard Airport air temperature record and other sites on the west coast are likely incorporated into the reanalysis, the quality of the reanalysis in the pre-satellite era is possibly even lower in remote areas with no observations. Due to the sparsity of available data in the overlap period 1979–2002, we have not evaluated the performance of Sval_Imp for variables other than air temperature. However, Østby et al. (2017) have evaluated the effect of this dataset discontinuity by simulating glacier mass balance using both ERA-40 and ERA-Interim to investigate whether this discontinuity could be responsible for a notable drop in simulated mass balance around the year 1980. They found that the ERA-40-based simulation yields an about 0.13 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ higher mass balance than the one based on ERA-Interim, but ERA-40-based simulations still show a 0.2 m drop of mass balance between 1970 and 1990, larger than that caused by the dataset discontinuity. This suggests that the change in mass balance regime was not caused by the heterogeneity of our composite forcing. Nevertheless, we cannot rule out the possibility that this change was caused by the discontinuity inherent in both reanalyses due to the availability of satellite observations after 1979. The annual observed air temperature trend for the period 1957–2013 at Svalbard Airport is 0.70±0.22 ${}^{\circ }\mathrm{C}{\mathrm{decade}}^{-\mathrm{1}}$, while the downscaled ERA data have an insignificantly lower warming trend of 0.67±0.19 ${}^{\circ }\mathrm{C}{\mathrm{decade}}^{-\mathrm{1}}$ at Svalbard Airport.

3.3 Radiation

To evaluate the quality of Sval_Imp radiation components, we use available records from two stations, roughly 300 km apart from each other. The record from Ny-Ålesund is from a daily serviced Baseline Surface Radiation Network station (Maturilli et al., 2015), whereas the Austfonna measurements are collected by an autonomously recording weather station (Schuler et al., 2014).

In Ny-Ålesund the model largely reproduces observations both for short and longwave radiation (Figs. 5 and 4, Table 1). During winter, downwelling longwave radiation is slightly underestimated, while there is no bias during summer. Since there is no air temperature bias in Ny-Ålesund during winter, the underestimation of longwave radiation is indicative of a too thin cloud cover in the reanalysis. In general, the representation of clouds are among the major issues of the reanalysis (Aas et al., 2016). Downwelling shortwave radiation is overestimated by 7 W m⁻² over the summer season in Ny-Ålesund. There is a much better agreement with radiation observations in Ny-Ålesund than on Etonbreen (Fig. 1), in northeastern Svalbard. This is to be expected, since radio soundings and other observation data from Ny-Ålesund are assimilated into ERA-Interim. Therefore, cloud cover at Ny-Ålesund is much better represented by the reanalysis than at Austfonna. On Etonbreen during summer, Sval_Imp underestimates downwelling shortwave radiation by 40 W m⁻², while downwelling longwave radiation is overestimated by 12 W m⁻², both indicative of an atmosphere that is too thick or too many clouds in the reanalysis. However, these biases could also be partly explained by measurement uncertainty caused by rime on the sensor or by sensor tilt. The latter issue is caused by the fact that the ice foundation of an autonomous weather station may melt and deform, causing tilt and thereby large errors, especially at high solar zenith angles (Bogren et al., 2016).

Table 1Meteorological stations used for validation of the downscaled reanalysis and their elevation (second line). N indicates the number of daily averages used in the validation, except for precipitation for which monthly sums have been evaluated. Seasonal biases in meteorological variables (downscaled minus observational averages) at all sites are averaged over each site's observation period. Shown are air temperature, T (K); relative humidity, RH (%); wind speed, WS (ms⁻¹); shortwave radiation (SW) and longwave radiation (LW) (both in W m⁻²); and precipitation, P (mm). The column headings S and W denote summer (June–August) and winter (September–May), respectively. Positive numbers indicate that the model results are larger than the observations. The second row at each site, shown in italics, is the bias between the raw ERA data and the observations.

^a Number of months. ^b Station located on a glacier.

Download Print Version | Download XLSX

3.4 Relative humidity and wind speed

For relative humidity the reanalysis represents the seasonality well, and in late summer both the humidity and the biases are of the largest magnitude. At the two coastal stations at Hopen and Rijpfjorden, the downscaled reanalysis is too dry, whereas it is too humid at the two higher-elevation stations. The coarse land mask of the reanalysis and the poor representation of sea ice are most likely the main causes for these biases. Wind speeds are reproduced reasonably well, including the seasonal cycle (Table 1). Biases are within ±1.5 m s⁻¹ with no clear seasonal trend. It is likely that the biases are caused by site-specific effects, such as deceleration of air flow in the lee of a topographic obstacle or acceleration due to channelizing through valleys.

Table 2Overview file structure of the Sval_Imp_v1 dataset.

Download Print Version | Download XLSX