Sval_Imp: A gridded forcing dataset for climate change impact research on Svalbard

. We present Sval_Imp, a high resolution gridded dataset designed for forcing models of terrestrial surface processes on Svalbard. The dataset is deﬁned on a 1km grid covering the archipelago of Svalbard, located in the Norwegian Arctic (74-82 ◦ N). Using a hybrid methodology combining multi-dimensional interpolation with simple dynamical modelling, the atmospheric reanalyses ERA-40 and ERA-interim by the European Centre for Medium-Range Weather Forecasting have been downscaled to cover the period 1957-2017 at steps of 6h. The dataset is publicly available from a data repository. In this paper, 5 we describe the methodology used to construct the dataset, present the organization of the data in the repository and discuss the performance of the downscaling procedure. In doing so, the dataset is compared to a wealth of data available from operational as well as project-based measurements. The quality of the downscaled dataset is found to vary in space and time, but generally represents an improvement compared to unscaled values, especially for precipitation. Whereas operational records are biased to low-elevations around the fringes, we stress the hitherto under-used potential of project-based measurements at higher elevation 10 and in the interior of the archipelago for evaluating atmospheric models. For instance, records of snow accumulation on large ice masses may represent measures of seasonally-integrated precipitation in regions sensitive to the downscaling procedure, thus providing added value. adapt the work, as long as you attribute the work to its origin and abide by the other license terms.


Introduction
The non-linearity of many surface processes poses challenges on appropriateness of atmospheric forcing for impact studies in terms of accuracy and precision (e.g., Liston and Elder, 2006). Especially in mountainous areas, the variability of surface 20 systems is typically governed by spatial scales not resolved in regional climate models and adjustments have to be made to overcome this (e.g., Fiddes and Gruber, 2014). A variety of methods has been developed for this purpose, differing in terms of data requirements and computational cost. While empirical-statistical scaling requires reference data for training and assumes a temporal robustness of the employed statistical relations (e.g., Ehret et al., 2012;Maraun, 2013), dynamic downscaling by means of high-resolution atmospheric modelling has high computational costs (e.g., Gutmann et al., 2016). 25 In this paper, we present Sval_Imp (Schuler, 2018), a high resolution gridded dataset obtained using a hybrid methodology combining multi-dimensional interpolation with simple dynamical modelling. The dataset is defined on a 1km grid covering the archipelago of Svalbard, located in the Norwegian Arctic (74-82 • N, 10-35 • E). The atmospheric reanalyses ERA-40 and ERA-interim by the European Centre for Medium-Range Weather Forecasting have been downscaled to cover the period 1957-2017 at 6 hours temporal resolution. The dataset comprises the near-surface variables required to compute the surface energy 30 balance, namely air temperature, precipitation, relative humidity, wind speed, and downwelling components of shortwave (solar) and longwave (thermal) radiation. Sval_Imp is publicly available from the Norwegian Research Data Archive NIRD, a data repository (https://doi.org/10.11582/2018.00006). In the following, we describe the methodology used to derive the dataset, present the organization of the data in the repository and discuss the performance of the downscaling procedure. For the latter, the dataset is compared to a wealth of data available from long-term operational as well as short-term scientific 35 records of meteorological and glaciological measurements. Operational records are biased to low-elevations around the fringes of the archipelago. Therefore, we stress the hitherto under-used potential of project-based measurements in the interior, highelevation regions for evaluating atmospheric models. For instance, records of snow accumulation on large ice masses may represent measures of seasonally-integrated precipitation in regions sensitive to the downscaling procedure, thus providing added-value. Sval_Imp has been employed entirely or in parts by a range of projects for forcing process models of the surface 40 energy and mass balances of glaciers (Østby et al., 2017), precipitation patterns and meltwater production in the Kongsfjord area (Pramanik et al., 2018), for assimilation of remotely-sensed snow cover using a snow distribution model (Aalstad et al., 2018) as well as for assessing growing conditions for fungi (Botnen et al., in prep.). Further, the dataset has been used to assess changes and trends in climate conditions of Svalbard (Hanssen-Bauer et al., 2019).

45
To generate fields of near-surface air temperature, precipitation, relative humidity, wind and downwelling shortwave and longwave radiation, we have downscaled the ERA-40 and ERA-interim reanalyses of the European Centre for Medium-Range Weather Forecasts (Uppala et al., 2005;Dee et al., 2011). The reanalysis data is provided at 6-hour intervals and has been retrieved on a 0.75 • × 0.75 • spatial grid covering the periods 1957-2002(ERA-40) and 1979-2017. These data have been downscaled to a 1km grid covering the region of interest ( Fig. 1) using the scheme described in the following 50 to produce Sval_Imp, a high-resolution, 6-hour dataset of precipitation, temperature, relative humidity, windspeed as well as downwelling shortwave and longwave radiation fluxes. of meteorological stations are represented by red diamonds, the numbers refer to the station names: 1 -Etonbreen, 2 -Janssonhaugen, 3 -Gruvefjellet, 4 -Kapp Heuglin, 5 -Rijpfjorden, 6 -Svalbard Airport, 7 -Isfjord Radio, 8 -Verlegenhuken, 9 -Hornsund, 10 -Kvitøya, 11 -Holtedahlfonna, 12 -Ny-Ålesund, 13 -Hopen. The black crosses indicate the grid points of the ERA reanalyses, and the grey countourlines indicate the topography of Svalbard in the ERA reanalyses at 0, 100, 200 and 300 m asl.

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 55 of the actual topography in the large-scale model used for the reanalysis (Fig. 1), leading to an underestimate 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 asl.
In contrast, the highest elevation in our gridded topography map is 1600 m asl. The roughness of the actual topography that gave rise to the name of the main island 'Spitsbergen' (sharp tops) is not represented by the smoothed topography used for the 60 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 which 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 70 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 Sections2.3 and 2.4.

75
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 80 of mountain wave to derive a transfer function that, for given wind conditions, relates the orographically enhanced precipitation to terrain topography.
By spectral decomposition and algebraic manipulation, Smith and Barstad (2004) derived the following transfer function .
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 σ(k, l) = U k + V l (U and V being the east and north components of the wind vector), and the conversion and fallout time scales τ c and τ f , respectively. In Equation (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.

90
Precipitation rates are obtained by retransformingP and adding it to the background precipitation P ∞ , that 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 Equation 1 to the largescale topography h ERA and removing the result from the ERA precipitation P ∞ = P ERA − P oro (h ERA ) (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.

100
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 of N m , τ f and τ c . Values of N m are calculated following where g is gravitational acceleration and T is vertically averaged air temperature weighted by the moisture content at several 105 pressure levels (Jarosch et al., 2012). Environmental lapse-rates Γ e are derived from air temperature at 700 hPa 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 −1 , representing an atmosphere with no stratification, and 0.01 s −1 representing a 110 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 time scales τ c and τ f are assumed equal and τ = τ c = τ 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 This setup is then applied to each 6h timestep and the resulting timeslices are progressively added to the record.

Air temperature
The downscaling for near-surface air temperature at 2 m level (T 2) closely follows the TopoSCALE procedure Fiddes and Gruber (2014), thereby we assume that the vertical structure of the free atmosphere determines the distribution of T 2 with 125 terrain elevation. For each 6h-time step, T 2 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 ERA field where T 2 > 273 K, instead of interpolating a 130 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/ ice surface the temperature of which is capped at 273 K. To avoid this effect, we assume that where T 2 > 273 K, the T 2 of the reanalysis is consistent with a melting surface and hence more realistic than a free-atmosphere interpolation.

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 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 140 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 cal-145 culated 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 emmissivity. This emissivity is estimated by accounting for a clear-sky component that depends on humidity and air temparature 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 150 multiplication with the sky-view factor.

Relative humidity and windspeed
Similar to our assumption about T 2 over a melting surface, we suggest that windspeed 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 (2017) for details. Furthermore, we present additional evaluation of the precipitation using snow measurements. Typically, 160 meteorological records are available as daily mean values and 6-hourly, downscaled variables have been temporally aggregated to match the time step of measurements. Precipitation in reanalysis is not constrained by data assimilation, giving rise to uncertainty in timing as well as amount. Our downscaling aims for reducing the bias in precipitation amount but does not treat the timing. For a performance evaluation, we therefore use monthly precipitation sums which are regarded as robust against timing mismatches, whereas a higher temporal resolution could penalize a method that otherwise is successful in reducing the 165 underestimation of the reanalysis.

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-  who evaluated the performance at 6 weather stations for several 30-year reference periods (1961-1990, 1971-2000, and 1988-2017). These biases are partly caused by too low precipitation measurements, 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 large undercatch for solid precipitation, however, this corrections 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 NE 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 Oct-May), they provide useful information about spatial precipitation patterns, in areas and elevations not covered by the operational meteorological stations. This approach builds 180 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 m w.e year −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 SE-NW asymmetry (Taurisano et al., 2007;185 Schuler et al., 2007;Dunse et al., 2009) caused by orographic enhancement of precipitation coming from the SE 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 190 shows almost no spatial variation and considerably underestimates observed values (Fig. 2). Østby et al. (2017) found that the winter mass balance of Hansbreen was not well reproduced both in terms of spatial pattern as well as accumulation amount. Aas et al. (2016) similarly reported the lowest performance for Hansbreen although they had 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

195
Hansbreen where the spatial distribution of snow is caused by wind redistribution rather than by the spatial precipitation pattern (Grabiec et al., 2006).

Air temperature
Downscaled air temparatures (Fig. 3) are compared to observations at the meteorological stations listed in Table 1 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 temparatures 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 • C and 1.5 • C at Svalbard Airport for ERA-40 and ERA-interim, respectively. We attribute the lower performance 210 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 at 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 sparsity of available data in the overlap period 1979-2002, we have not evaluated the performance of Sval_Imp for other variables than air temperature. However, Østby et al. (2017) have evaluated the effect of this dataset discontinuity by simulating glacier mass in simulated mass balance around year 1980. They found that the ERA-40 based simulation yields an about 0.13 m w.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 data set 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 220 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 • C decade −1 , while the downscaled ERA data has an insignificantly lower warming trend of 0.67±0.19 • C decade −1 at Svalbard Airport.

Radiation
To evaluate the quality of Sval_Imp radiation components, we use available records from two stations, roughly 300 km apart 225 from each other. The record from Ny-Ålesund is from a daily serviced Baseline Surface Radiation Network station (Maturilli et al., 2013) 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, Tab. 1).
During winter, downwelling longwave radiation is slightly underestimated, while there is no bias during summer. Since there 230 is no air temparature 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 −2 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.

235
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 −2 while downwelling longwave radiation is overestimated by 12 W m −2 , both indicative of a too thick atmosphere 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 stations may melt and deform causing tilt and 240 thereby large errors, especially at high solar zenith angles (Bogren et al., 2016).

Relative humidity and windspeed
For relative humidity the reanalysis represents the seasonality well, and in late summer both the humidity and the biases are of 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 245 most likely the main causes for these biases.
Wind speeds are reproduced reasonably well, including the seasonal cycle (Tab. 1). Biases are within ±1.5 m s −1 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.    The downloadable dataset comprises individual files for each of the variables precipitation, air temparature, relative humidity, org/index.php/Attribute_Convention_for_Data_Discovery_1-3). Table 2 gives an overview over the number of files and their sizes for the different epochs (ERA-40, ERA-interim) and 265 variables contained in the dataset. The naming convention for the individual files is <EPOCH>_<VAR>_<YYYYMM>.nc, where EPOCH is either "ERA40" or "ERAi", VAR is an abbreviation of the variable of interest (one of "precip", "temp", "RH", "windspeed", "SWi" or "LWi") and YYYYMM identifies year and month, for instance ERAi_temp_200410.nc is the name of the file containing air temparature from ERA-interim for October 2004.

270
We present a gridded dataset of near-surface, meteorological variables at 1km resolution covering the Svalbard archipelago.
The set of variables enables application of energy balance models and comes at a time steps of 6 h. The high-resolution grids are derived from coarse-scale reanalyses ERA-40 for the period 1957-2002 and from ERA-interim for 1979-2017. We describe