Merits of novel high-resolution estimates and existing long-term estimates of humidity and incident radiation in a complex domain

To provide better and more robust estimates of evaporation and snowmelt in a changing climate, hydrological and ecological modeling practices are shifting towards solving the surface energy balance. In addition to precipitation and near-surface temperature (T₂), which often are available at high resolution from national providers, high-quality estimates of 2 m humidity and surface incident shortwave (SW↓) and longwave (LW↓) radiation are also required. Novel, gridded estimates of humidity and incident radiation are constructed using a methodology similar to that used in the development of the WATCH forcing data; however, national 1 km×1 km gridded, observation-based T₂ data are consulted in the downscaling rather than the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ Climatic Research Unit (CRU) T₂ data. The novel data set, HySN, covering 1979 to 2017, is archived in Zenodo (https://doi.org/10.5281/zenodo.1970170). The HySN estimates, existing estimates from reanalysis data, post-processed reanalysis data, and Variable Infiltration Capacity (VIC) type forcing data are compared with observations from the Norwegian mainland from 1982 through 1999. Humidity measurements from 84 stations are used, and, by employing quality control routines and including agricultural stations, SW↓ observations from 10 stations are made available. Meanwhile, only two stations have observations of LW↓. Vertical gradients, differences when compared at common altitudes, daily correlations, sensitivities to air mass type, and, where possible, trends and geographical gradients in seasonal means are assessed. At individual stations, differences in seasonal means from the observations are as large as 7 ^∘C for dew point temperature, 62 W m⁻² for SW↓, and 24 W m⁻² for LW↓. Most models overestimate SW↓ and underestimate LW↓. Horizontal resolution is not a predictor of the model's efficiency. Daily correlation is better captured in the products based on newer reanalysis data. Certain model estimates show different dependencies on geographical features, diverging trends, or a different sensitivity to air mass type than the observations.

1 Introduction

Geophysical modeling is advancing, and more and more hydrological, ecological, and land surface models (from here on referred to as land models) are now estimating the surface energy balance (Mueller et al., 2013). Shortwave radiation is the exogenous energy provider to Earth. At middle and higher latitudes, surface downward longwave radiation is an equally important radiative driver at the surface. Estimating the surface energy balance provides a sensible heat flux, as well as a latent heat flux, which in turn can be converted to evaporation or snowmelt, key variables for estimating the surface water balance.

Recent studies have shown the added value of using more forcing data than only precipitation and temperature when modeling evaporation (Milly and Dunne, 2011; Lofgren et al., 2011; Haddeland et al., 2012; Pierce et al., 2013; Stagge et al., 2014) or snow cover (Raleigh et al., 2016; Harpold et al., 2017). High-quality and robust diagnoses, forecasts, and projections of evaporation- and snowmelt-related processes are essential for flood and hydropower management. Further, gridded data sets of high quality are needed to statistically bias-correct or downscale future climate scenarios (Abatzoglou, 2013), to spin up land surface models (e.g., Rodell et al., 2005; Koster et al., 2004; Kristiansen et al., 2012), and to assist model development.

Considerable effort is used to improve process description in environmental models and compare the results of different models. When land models are run without coupling to an atmospheric model, i.e., in offline or stand-alone mode, meteorological near-surface variables, commonly referred to as forcing data, are required. In practice, different communities use different forcing data estimates, such as the more empirically based estimates from the MTCLIM algorithms (Bristow and Campbell, 1984; Thornton and Running, 1999; Bohn et al., 2013), estimates from numerical weather and climate models, or a combination of the two (see, e.g., Mizukami et al., 2016). The different approaches used make it more difficult to compare the output across land models and have resulted in dedicated projects where various models are run with similar forcing data in, for example, the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP; Warszawski et al., 2014).

The Norwegian operational hydrological models have historically been calibrated and adapted to use high-resolution, gridded 2 m temperature and precipitation data as forcing data. For this purpose high-resolution (1 km×1 km) data sets, which cover the period from 1957 until the present (2019) with a daily resolution have been developed, namely the SeNorge data (Mohr, 2008; Tveito and Førland, 1999; Lussana et al., 2018a, b). In recent times a long-term high-resolution, quantile-mapping-based gridded data set of near-surface wind speed has also been developed at the Norwegian Meteorological Institute (MET Norway, available from http://thredds.met.no/thredds/catalog/metusers/klinogrid/KliNoGrid_16.12/FFMRR-Nor/catalog.html, last access: 10 June 2019). Gridded, observation-based, high-resolution data sets for humidity and incident radiation are, however, lacking.

Previous studies have compared and validated gridded estimates of humidity or incident radiation globally (Bohn et al., 2013; Schmied et al., 2016; Weedon et al., 2011), for regions in the US (Slater, 2016; Mizukami et al., 2014; Pierce et al., 2013; Lapo et al., 2017), coastal Brazil (Almeida and Landsberg, 2003), France (Szczypta et al., 2011), and the pan-Arctic region (Shi et al., 2010). No previous studies have, as far as we know, compared and assessed the quality of high-resolution empirically based and reanalysis-based estimates of humidity and incident radiation for regions within Europe (let alone for Norway specifically). This results in an additional and unnecessary source of uncertainty for land modeling in Norway.

Norway's complex topography and coastline may suggest that high-resolution data sets would perform better. When land models are run the near-surface temperature from the forcing data is usually adjusted to sea level and then to the land model's fine-resolution grid cell elevation that uses a standard atmospheric lapse rate to account for the difference in terrain height in the forcing data model and the land model to be run. A standard atmospheric temperature lapse rate may be unreasonable in winter, at high latitudes (Kotlarski et al., 2010; Brinckmann et al., 2016), and in complex terrain (Mizukami et al., 2014), and a much lower resolution in the forcing data grid compared to the land model may increase these error components. Further, if humidity is not also adjusted for inconsistencies between temperature and humidity will likely result in an unrealistic relative humidity (Haddeland et al., 2006; Weedon et al., 2011). Incident radiation is influenced by variables showing a strong vertical dependence like near-surface temperature and humidity, cloud cover (e.g., Marty, 2000), and local variations in surface components like vegetation and snow cover (Erlandsen et al., 2017; Rydsaa et al., 2017), and thus may likely benefit from vertical adjustment.

While the spatial correlation may improve in a data set with a high spatial resolution, Decker et al. (2012) highlight the need to address temporal correlation on timescales shorter than monthly scales in data constructed from reanalyses for the purposes of forcing land surface models. A high horizontal resolution may lead to a better representation of the average state of a variable but not necessarily to an improved description of the concurrent temporal evolution of the forcing variables on shorter timescales, e.g., during the rapid passage of a low-pressure system with multiple distinct air mass characteristics and precipitation types.

This study addresses the aforementioned sources of uncertainty concerning commonly used estimates of humidity, either in the form of vapor pressure (VP) or converted to dew point temperature (T_d), incident longwave radiation (LW↓), and incident (global) shortwave radiation (SW↓) available for long-term land surface modeling in the region through the following processes.

• The construction of an original data set, HySN (to explore the benefit of utilizing a 1 km×1 km national data set of 2 m temperature in the post-processing reanalysis data);

• The gathering of global long-term gridded data sets of humidity and incident radiation from two reanalysis data sets, two post-processed reanalysis data sets, and two versions of empirically based estimates compiled for continental Norway;

• The aggregation of available observations of humidity and incident radiation between 1982 and 1999 from a variety of providers, and where necessary, implementing quality control routines;

• The construction of multiple linear regression models to provide vertical gradients in both the observations and the model estimates, so that the variables may be adjusted to a similar altitude before their differences are assessed;

• The correlation of model estimates with observations on a daily timescale is explored by compiling anomaly correlation coefficients;

• The comparison of the model estimates' cumulative distributions; their sensitivity to weather types, continentality, and latitude; and their decadal trend to the observational data.

Additionally, two hypotheses are investigated. ℋ_a – there are vertical gradients in near-surface humidity and incident radiation in our domain. ℋ_b – the added value of the high horizontal resolution of the more empirically based estimates outweighs the added value of relying on estimates from coarser-resolution numerical weather prediction reanalyses.

The data sets considered are two global reanalysis data sets, the NASA Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA2) (Bosilovich et al., 2015, 2017), and ECMWF's ERA-Interim (Dee et al., 2011), two products based on reanalysis data post-processed using higher-resolution gridded observational data, the Princeton Global Meteorological forcing data set, version 2 (PGMFDv2) (Sheffield et al., 2006), the WATCH forcing data methodology applied to ERA-Interim (WFDEI) (Weedon et al., 2014), and two versions of high-resolution empirically based estimates from the preprocessor of the Variable Infiltration Capacity (VIC) model, a macroscale hydrological model (Liang et al., 1994) largely based on the MTCLIM algorithms. Finally, a novel data set, the HySN data set, is compiled for the current study and evaluated. HySN is compiled by employing a similar method as was used in the development of PGMFDv2 and WFDEI; however, in this case ERA-Interim near-surface humidity and incident radiation are post-processed using a national, high-resolution, gridded 2 m temperature data set, SeNorge2 (Lussana et al., 2018b).

2 The gridded humidity and radiation estimates considered

Long-term data sets that are freely available, which can be used to drive hydrological, ecological, and land surface models for the Norwegian domain, include the newer reanalyses: MERRA2 and ERA-Interim. Due to computational constraints, currently available long-term global reanalysis data have horizontal resolutions ranging from $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ to $\mathrm{0.5}{}^{\circ }\times \mathrm{0.66}{}^{\circ }$. The MERRA2 reanalysis has a resolution of 0.5^∘ latitude $\times \mathrm{0.66}{}^{\circ }$ longitude. Around Oslo, Norway, this corresponds to a grid cell height and length of about 56 km×42 km. The reanalysis data sets are based on global circulation models ingesting large amounts of observational data by making use of complex assimilation techniques. However, substantial biases may still occur in reanalysis data. Heikkilä et al. (2011) found a mean error of +42.9 % in precipitation intensity in ERA-40 over Norway between 1961 and 1990. Bromwich et al. (2016) found a negative bias in ERA-Interim surface LW↓ radiation and precipitation between November 2007 and December 2008 across middle and high latitudes in the Northern Hemisphere.

The coarse resolution of reanalysis data and the knowledge of biases that may be present in them has spurred the development of post-processed reanalysis data sets. The PGMFDv2 and WFDEI are data sets consisting of variables relevant for forcing land surface models. The relevant variables are extracted from reanalysis data and post-processed and downscaled with gridded observational data. Both data sets are global and have horizontal resolutions of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$. PGMFDv2 and WFDEI both adjust reanalysis estimates of humidity and LW↓ with the gridded, $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ Climatic Research Unit (CRU) T₂ following the methods described in the development of NLDAS (Cosgrove, 2003). Taking a note from these methods, a novel high-resolution product is developed and validated in the current study: Hybrid SeNorge, abbreviated as HySN. HySN is constructed by post-processing ERA-Interim humidity and radiances in a similar manner to PGMFDv2 and WFDEI but utilizing a national data set, the 1 km×1 km SeNorge2 T₂, rather than the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ CRU T₂.

Another source of near-surface humidity and incident radiation estimates are the MTCLIM algorithms, which combine first principles from atmospheric physics with empirical extrapolation logic. Precipitation and temperature, variables that often are available from a dense network of surface observation stations, are used to estimate shortwave radiation and humidity. Versions of the MTCLIM routines are used to provide forcing data for a large number of hydrological and ecological models; it has, for example, recently been made available for the Mesoscale Hydrological model (MHm v5.9, https://doi.org/10.5281/zenodo.1069202). The variables estimated from MTCLIM are often utilized for impact studies, e.g., the impacts of climate change and forest management on ecosystem services in Europe (Bugmann et al., 2017). The algorithms have also been used to generate several gridded data set products of humidity and radiation for the US (e.g., Livneh et al., 2013) and are used to provide climate change projections of humidity and radiation for the US (Bureau of Reclamation, 2013). However, several recent studies have found regionally inconsistent biases in the MTCLIM algorithms (Shi et al., 2010; Bohn et al., 2013; Pierce et al., 2013; Slater, 2016; Mizukami et al., 2014).

The orography and land masks of the models are presented in Fig. 1. Compared to ERA-Interim orography, the SeNorge grid elevation is on average higher (mean: 37 m, median: 13 m; see the red areas in Fig. 1). The difference in maximum elevation is more than 1000 m. Meanwhile, near the coast and in inland areas the ERA-Interim orography is predominantly higher (see the blue areas in Fig. 1). The data sets are summarized in Table 1. Further details considering ERA-Interim, MERRA2, WFDEI, PGMFDv2, WFDEI, HySN, and two data sets from the VIC model's preprocessor, largely based on the MTCLIM algorithms, are presented in the following.

Cosgrove (2003)Cosgrove (2003)Cosgrove (2003)Thornton and Running (1999)

Table 1 The following data sets provide estimates of humidity, LW↓, and SW↓, which are then evaluated. Precipitation is denoted as P, 2 m temperature as T₂. The global data sets are retrieved from online repositories, while the data sets with regional coverage are compiled locally based on the stated input data, the VFD data sets using the VIC model's preprocessor, and the HySN data set based on the methods outlined in the current study. Additional references are given elsewhere in the text.

2.1 ERA-Interim

The ERA-Interim (Dee et al., 2011) is a reanalysis data set developed by ECMWF, covering the time period from 1979 until the present. It is based on a 2006 release of the ECMWF operational model system (IFS Cy31r2) and has a horizontal resolution of about 80 km. It includes a four-dimensional variational analysis (4D-Var). Surface observations are ingested by optimal interpolation. The variables evaluated in this study are daily means of 2 m temperature and dew point temperature from analysis times (00:00, 06:00, 12:00, and 18:00 UTC) and LW↓ and SW↓ taken between +12 and +24 h into the forecast to allow for spin-up (see, e.g., Weedon et al., 2014).

2.2 Modern-Era Retrospective Analysis for Research and Applications 2 (MERRA2)

MERRA2 is an atmospheric reanalysis data set developed by NASA, available from 1980 until the present, with a horizontal resolution of 0.5${}^{\circ }\times \mathrm{0.625}$^∘(Bosilovich et al., 2015). Mass conservation constraints are imposed so that assimilated observations do not cause large violations of the global water balance. In MERRA2 land surface observations are not assimilated. The data variables used in this study are model orography (Mer, 2015a), pressure, and humidity from atmospheric single level diagnostic (Mer, 2015c), LW↓, and SW↓ (Mer, 2015b).

2.3 Princeton's global meteorological forcing data set version 2 (PGMFDv2)

PGMFDv2 is an updated version of the 0.5${}^{\circ }\times {\mathrm{0.5}}^{\circ }$ 60-year Princeton global meteorological forcing data set (Sheffield et al., 2006). The updates are described in Schmied et al. (2016). The humidity, LW↓, and SW↓ estimates are based on the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP-NCAR) reanalysis but post-processed to comply with the gridded, observation-based time series of precipitation, temperature, and cloud cover, with a horizontal resolution of 0.5${}^{\circ }\times {\mathrm{0.5}}^{\circ }$ from the Climatic Research Unit (CRU TS 3.2.1) and satellite estimates of LW↓ and SW↓.

2.4 The WATCH forcing data methodology applied to ERA-Interim (WFDEI)

The application of the WATCH forcing data methodology to ERA-Interim reanalysis data, WFDEI, is described in Weedon et al. (2014). The data are available from 1979 to the present and have a horizontal resolution of 0.5${}^{\circ }\times {\mathrm{0.5}}^{\circ }$. The humidity, LW↓, and SW↓ estimates are based on ERA-Interim data, post-processed to comply with the global gridded, observation-based time series of 2 m temperature, cloud cover, and interannual aerosol loading from CRU TS, using CRU TS 3.2.1 prior to 2009 similar to PGMFDv2.

2.5 Hybrid SeNorge ERA-Interim, HySN (H)

As part of this study, additional estimates of humidity and LW↓ are derived, using methods based on Cosgrove (2003), adjusting ERA-Interim humidity and LW↓ to comply with the newly developed, 1 km×1 km SeNorge2 T₂ data set (Lussana et al., 2018b). Further, the ERA-Interim SW↓ estimates are adjusted based on the previously adjusted humidity estimates and the 1 km ×1 km orography.

ERA-Interim humidity and longwave radiation are vertically adjusted on a daily basis by consulting the daily SeNorge T₂. The method differs from that used in the construction of the WFDEI and Princeton forcing data, where the reanalysis T₂ is adjusted to sea level and then to the CRU grid elevation using a constant lapse rate, before adjusting it on a monthly basis to fit the monthly mean CRU T₂. The vertical adjustment of humidity makes use of the common approximation that relative humidity remains constant with height (see, e.g., Feld et al., 2013), making it easy to solve for a SeNorge2 compatible dew point temperature T_d based on ERA-Interim relative humidity (RH) and SeNorge2 T₂. Humidity is corrected to saturation if supersaturation occurs. Surface pressure is adjusted using an approximation of the hypsometric equation. The vertical adjustment of longwave radiation is done by scaling an empirical expression for clear-sky LW↓ to the SeNorge T₂ and the previously compiled vertically adjusted humidity estimate. No consistent approach is used in other forcing data sets when vertically adjusting SW↓. Given that SW↓ is very sensitive to near-surface humidity and that the Cosgrove (2003) method used above adjusts humidity, we chose to scale the ERA-Interim SW↓ to the ratio of estimated clear-sky transmissivity calculated using an empirical equation from Thornton and Running (1999), taking into account the difference in altitude and humidity in the two data sets. A clear-sky type correction approach is thus used to adjust both SW↓ and LW↓.

For consistency with SeNorge precipitation and T₂, the variables have a temporal resolution of a day, starting from 06:00 UTC. The data are currently compiled for the time period 1979–2017 and cover the same domain as the SeNorge2 grid. The data compilation is described in detail in the Supplement. The HySN data product is freely available from Zenodo (https://doi.org/10.5281/zenodo.1970170), and the Python code to generate the data is available on GitHub (https://doi.org/10.5281/zenodo.1435555).

2.6 VIC type forcing data, VFDv1, and VFDv2

The humidity and radiation estimates referred to here as VIC type forcing data (VFD) (see, e.g., Bohn et al., 2013; Pierce et al., 2013) are products of the VIC model's preprocessor. The VIC model includes the option to generate gridded humidity and radiation from gridded daily precipitation and maximum and minimum temperature. The VIC model preprocessor includes a slightly modified version of the MTCLIM model and algorithms for estimating longwave radiation. The MTCLIM algorithms included in the VIC preprocessors to estimate humidity use a modified version of the Magnus formula with daily minimum temperature used as a proxy for T_d (Kimball et al., 1997). Shortwave radiation is estimated using the Thornton and Running algorithm (Thornton and Running, 1999). The variables are estimated simultaneously; i.e., the algorithms supply each other with information (Bohn et al., 2013).

Two versions of VFD are evaluated in this study. The first version, from here on called VFDv1, uses daily precipitation and mean temperature from SeNorge version 1.1 (Tveito and Førland, 1999; Mohr, 2008), supported by hourly temperature fields from a regional atmospheric reanalysis data set, NOrwegian ReAnalysis (NORA10, Reistad et al., 2011), with a resolution of about 11 km to compile maximum and minimum temperature using a method similar to Vormoor and Skaugen (2013). The VIC4.0.6 preprocessor is used with default options, i.e., a modified version of MTCLIM4.2, and the TVA clear-sky and Bras full-sky LW↓ algorithm (Bras, 1990).

The second version of VIC type forcing data, from here on called VFDv2, is based on slightly different input data, i.e., precipitation and mean, maximum, and minimum temperature from a newer version of the 1 km by 1 km SeNorge data, SeNorge2 (Lussana et al., 2018a, b). The VIC4.0.6 preprocessor is used with default options, i.e., a modified version of MTCLIMv4.3, and with LW↓ estimates based on the Prata (1996) clear-sky algorithm combined with the Deardorff (1978) full-sky algorithm (for further references see Bohn et al., 2013).

3 Study area

Norway is located in the receiving end of the westerlies that pass over the North Atlantic. This, combined with a long coast lined with mountains provides Norway with 1500 mm of precipitation a year, with distinct regional differences in precipitation amounts received. Although almost 40 % of Norway is covered by forest, evaporation from the land surface is estimated to be less than a fourth of the received precipitation (Hanssen-Bauer et al., 2009). Most of Norway will normally have snow cover in the winter season, with the length of the snow season varying from a few days to 300 d a year (dependent on latitude, elevation, and distance from the coast). Mean temperature (1971–2000) is 1.3 ^∘C, and varies from 7 ^∘C near the coast in southern Norway to −4 ^∘C in the mountains. Between 1976 and 2014 T₂ increased by half a degree Celsius per decade (Hanssen-Bauer et al., 2017).

4 Surface observations

The model estimates and observational data are compared between 1982 through 1999. The observational data include humidity measurements from 84 sites, SW↓ observed at 10 sites, and LW↓ observed at two sites. The comparison of the model estimates of incident radiation with stations data is only made possible by including observations from agricultural stations and applying quality control routines. The observations are gathered from the University of Bergen (UiB, SW↓ and LW↓ measurements), and from the Norwegian Meteorological Institute's repository for observational data, which also includes measurements from agricultural stations conducted by the Norwegian Institute of Bioeconomy Research (NIBIO). The locations of the stations used in the comparison with the model estimates are shown in Fig. 1.

Figure 1 Panels (a)–(e) show the orography and land mask of MERRA2 (a), ERA-Interim (b), PGMFDv2 (c), WFDEI (d), and SeNorge (c), respectively, visualized on the SeNorge UTM33 grid with a green–brown color scheme. For reference, national borders and the coastline derived from a high-resolution data set are delineated in black. The locations of the 84 VP stations used in the model comparison are denoted with red crosses in (f). The locations of SW↓ and LW↓ stations are marked in (g) with red and orange markers, respectively. Note that the southernmost LW↓ station also measures SW↓. The last map (h) displays the difference in meters between the SeNorge and ERA-Interim orography in common land areas. Higher elevations in SeNorge are indicated with red, while blue indicates higher elevations in ERA-Interim.

4.1 Humidity

Humidity observations from 84 stations are included in the study (see Fig. 1). A minimum of 5 years of daily data was necessary for the station data to be included; however, most stations have the complete station record (18 years) available. The observations, once converted to vapor pressure, have an uncertainty ranging from around 5 % at 20 ^∘C to 6 % at −20 ^∘C (Gabriel Kielland, personal communication, MET Norway, 2019). The latitude, longitude, altitude, distance to the ocean, and the start and end date of the time series are given for each station in the Supplement (Tables S1 and S2).

4.2 Shortwave radiation

The location of the 10 stations included in the evaluation of modeled SW↓ is displayed in Fig. 1 and covers a latitudinal range from 58.8 to 69.7^∘ N. Seven of the stations are agricultural stations managed by NIBIO. The latitude, longitude, altitude, distance to the ocean, the start and end date of the time series, and the percentage of flagged data are given for each station in the Supplement (Table S3). The number of days of data used in the validation varies from 5.6 years for Gjengedal to more than 17 years for Bergen.

Most stations measure global radiation with a Kipp & Zonen CM11 thermoelectric pyranometer. The estimated uncertainty of hourly and daily totals of CM11 may be as low as 3 % in optimal conditions (Grini, 2015). Daily global SW↓ measurements from Bergen station (UiB) are included in the World Radiation Data Centre (WRDC) and are quality-controlled by the data provider (UiB, A. Olseth). The daily estimates have an uncertainty of 3.5 % (Parding et al., 2016). Measurement errors and uncertainty may depend on sensor calibration, placement (e.g., sky-view), the temporal resolution of the measurements, cleaning of the pyranometer, and local weather conditions.

For stations other than Bergen, quality control procedures follow the methodology suggested by Grini (2015) as outlined in Table 2. This procedure involves running rtmrun (Godøy, 2013); a Perl wrapper around Libradtran 1.7 (Mayer and Kylling, 2005); a library for radiative heat transfer to provide solar zenith angle (SZA), extraterrestrial (SW_E↓), and clear-sky (SW_CS↓) incident shortwave radiation for each station location; and the Python scripts developed in Grini (2015) to screen and flag the data based on automatic quality control tests. Measurements exceeding the upper and lower bounds given in Table 2 were flagged. Additionally, all station time series were visually inspected at hourly, daily, and monthly levels in order to flag erroneous data not captured by the automatic routines, with emphasis on data points marked as suspicious due to large hourly increments or very high or low variation in the ratio of observed to extraterrestrial radiation. Figure 2 shows the mean monthly values of SZA, SW_E↓, SW_CS↓, the raw measured values (SW_raw↓), and values passing the quality control routines (SW_QC↓) for Løken station.

Figure 2 Estimated and measured shortwave radiation (W m⁻²) at Løken station (61.1^∘ N, 9.1^∘ E) for the period January 1991 to December 1999. Mean monthly solar zenith angle (SZA), (global) shortwave incident radiation at the top of the atmosphere (SW_E↓), modeled clear-sky incident shortwave radiation SW_CS↓), and station measurements of incident shortwave radiation before (SW_raw↓) and after (SW_QC↓) quality control are shown.

Table 2 An overview of the automatic quality control tests, based on the relative values of the solar zenith angle (SZA), measured (SW_raw↓), extraterrestrial (SW_E↓), and clear-sky (SW_CS↓) incident global radiation. The table is adapted from Table 4.1.1 in Grini (2015).

In the calculation of daily means, values were flagged as erroneous and subsequently excluded from the validation if more than two hourly data points were flagged or missing during daytime. The number of discarded days varied from 4 % at Kise to 29 % at Tromsø.

4.3 Longwave radiation

Bergen station and Voll station (Trondheim), denoted with orange markers in Fig. 1g, have observations of incident longwave radiation available for the time period considered. The lack of LW↓ measurements is not an uncommon challenge (see, e.g., Carrer et al., 2012). The stations' latitude, longitude, distance to the coast, and the start and end date of the data used are listed in the Supplement (Table S3). Both measurement stations are located within 5 km of the coast (see Fig. 1). The measurements from Bergen are managed and quality-controlled by UiB. In the first part of the period they are from a Schulze radiation balance meter, while later in the period they are from an Eppley pyrgeometer. The sensors are placed on the roof of UiB. The observation station at Voll, Trondheim, was managed by MET Norway from March 1996, until it was shut down. The Trondheim measurements are from a Kipp & Zonen CG 1 pyrgeometer located on the ground. At both stations unshaded sensors were used, possibly leading to slight overestimation due to solar near-infrared radiation contamination (overestimations of 10 % on cloud-free days were found in de Oliveira et al., 2006; Meloni et al., 2012). The data were quality-controlled by visual inspection for spikes and jumps and by comparing the consistency between the two time series. The Stefan–Boltzmann blackbody longwave radiation was set as an upper limit of the measurements, using the air temperature from the station. If more than 2 h were missing or flagged during a day, observations from that day were omitted in the subsequent validation.

5 Evaluation methods

Daily estimated values are compared to station observations. The nearest model grid cell is selected from the data sets without interpolation, to avoid introducing spatial or temporal smoothing of the meteorological fields (see Hofstra et al., 2010; Gutmann et al., 2012). This study specifically looks into the altitudinal dependence of the humidity and surface incident radiation estimates, and as a starting point the estimates without adjustment to the observation stations' altitudes are used.

5.1 Vertical adjustment to station altitude

Prior to the comparison of the model estimates with the station observations, the observations and model estimates of VP and SW↓ were analyzed using multiple linear regression with geographical features as predictors, in order to find vertical gradients so that the model estimates could be adjusted to station altitude. The geographical predictors used were altitude (either the stations' altitude or, for the models, the altitude of the nearest-neighbor grid cell to the stations), latitude above 57^∘ N, and distance to the coast, which was calculated in Python using the Haversine distance from the station to the coastline, extracted from a coastline data set (Wessel and Smith, 1996) that is available via the Matplotlib Basemap Toolkit, implemented at a coarse resolution to not include large inland lakes. The limited number of SW↓ measurement stations and the varying temporal availability of high-quality observational data made an evaluation of the altitudinal dependence of the SW↓ more demanding. The SW↓ data were first converted to clearness index (CI), which describes the daily incident shortwave radiation fraction of the potential extraterrestrial radiation at the local position and time (SW$\downarrow /$ SW_E↓), and then daily data of over 1000 concurrent measurements from eight stations were used in the regression.

The model estimates were adjusted to station altitude by multiplying their grid cell values with the difference in altitude to the observation station and a vertical adjustment gradient. For each model the vertical adjustment gradients were computed as the mean of coefficients found for the model in question and those found for the observational data, linearly interpolated from a seasonal to a daily frequency. A similar regression model was constructed to find vertical gradients in LW↓ using a well-performing data set in lieu of observations due to the limited number of LW↓ observation stations.

5.2 Evaluation metrics

Seasonality and aggregated means are assessed by plotting the mean monthly station values for the observations and models. The differences between the model estimates and observations at individual stations are displayed in heat maps. For each variable a table is provided listing several metrics. The tables list the following variables for each model:

1. $\mathrm{\Delta }=\stackrel{\mathrm{‾}}{{\mathit{\mu }}_{\text{station,model}}-{\mathit{\mu }}_{\text{station,observation}}}$, the mean (μ) of the station differences;

2. $\left|\mathrm{\Delta }\right|=\stackrel{\mathrm{‾}}{\left|\right({\mathit{\mu }}_{\text{station,model}}-{\mathit{\mu }}_{\text{station,observation}}\left)\right|}$, the mean of the station absolute differences;

3. $|\mathit{\delta }{|}_{max}=max(|{\mathit{\mu }}_{\text{station,model}}-{\mathit{\mu }}_{\text{station,observation}}|)$, the largest absolute difference at any station;

4. $|{\mathit{\delta }}^{\mathrm{s}}{|}_{max}=max(|{\mathit{\mu }}_{\text{season,station,model}}-{\mathit{\mu }}_{\text{season,station,observation}}|)$, the largest absolute difference found at any station in any season.

Also listed are the mean daily anomaly correlation coefficient (ACC), i.e., the daily Pearson correlation coefficient of the time series where the observed day-of-year mean is subtracted and the number of stations where the cumulative distribution of daily mean estimates passes (p>0.001), and finally the Kolmogrov–Smirnov test of similarity with the cumulative distribution of the observations. The Kolmogrov–Smirnov test returns the probability that the underlying one-dimensional probability distributions are the same (H₀). The similarity between the models and observations on a daily frequency is visualized in Taylor plots, where the normalized standard deviation, the root-mean-square error, and the correlation coefficient of the de-seasonalized time series are displayed. The time series are de-seasonalized by subtracting the observed day-of-year climatology. The correlation coefficient thus corresponds to a non-centered version of the anomaly correlation coefficient (ACC).

5.3 Evaluation of geographical gradients

In order to see if the geographical dependencies of the model estimates of humidity and shortwave radiation differ significantly from those seen in the observational data, similar multiple linear regression models as those previously constructed to find vertical gradients are used. The predictors are the season, altitude (z), latitude above 57^∘ N, and distance to the coast (C). The regression is first performed separately for each model and for the observations. However, a second iteration of regression is performed for each of the models, where the input data are composed of the observational data and model data are appended together with the data sources indicated. The data source is then used as a categorical predictor allowed to interact with any of the model coefficients. Significance of the interaction term, e.g., between latitude and model source, will indicate that the model's latitudinal gradient is significantly different from the gradient seen in the observations.

5.4 Air mass type sensitivity

The differences between the model estimates and observations are also inspected for an air mass type dependence. Bower et al. (2007) found significant decreases in the frequencies of dry moderate and dry polar air mass types at Bergen (Flesland) between 1974 and 2000. If the precision of the model estimates is dependent on air mass type, the derived changes in the variables with time may be less robust if the frequencies of the air mass types also change with time. The spatial, synoptic air mass type classification has been constructed for 48 stations in Europe and 7 stations in Norway (Sola, Flesland, Fornebu, Ørlandet, Bodø, Tromsø, and Slettnes) by Bower et al. (2007), according to the methods developed in Sheridan (2002) and Kalkstein et al. (1996). The categorization is done by using sub-daily surface observations of temperature, dew point, wind, pressure, and cloud cover at individual stations (often airports). The synoptic weather-typing classifies the local air mass conditions into the categories DP (dry polar), DM (dry moderate), DT (dry tropical), MP (moist polar), MM (moist moderate), MT (moist tropical), and TR (transitional). The dry weather types are associated with clearer conditions, while the moist weather types are associated with clouds and higher humidity. TR days are defined by large shifts in the synoptic variables, i.e., days where the weather type is changing. The MT weather type is often found in the warm sectors of cyclones, while the MP and MM type may be found in the vicinity of a front or in air transported inland from a cool ocean.

5.5 Comparison of trends

The year 1985 is considered the start of the SW↓ brightening period in Europe after a period of SW↓ dimming due to aerosol emissions (see, e.g., Wild et al., 2005). For time series of a sufficient length and quality, the observational data and the model data are inspected for trends. Stations that have less than 10 % missing daily data between January 1985 and December 1999 are considered, and trends are calculated for each calendar month using the Mann–Kendall test and by calculating the Sen slope (Hirsch et al., 1982). The analysis is done using the R software and functions within the “trend” package (Pohlert, 2018). A total of 59 humidity stations meet the criteria and are grouped into five geographical regions, southwest (S_W, stations with  $<\mathrm{61.8}{}^{\circ }$ N, $<\mathrm{8}{}^{\circ }$ E), southeast (S_E, stations with  $<\mathrm{61.8}{}^{\circ }$ N,  $>\mathrm{8}{}^{\circ }$ E), central (C, $\mathrm{61.8}>{}^{\circ }$ N$<\mathrm{64}{}^{\circ }$ N), northwest (N_W, $>\mathrm{64}{}^{\circ }$ N, $<\mathrm{20.6}{}^{\circ }$ E), and northeast (N_E, >64^∘ N, >20.6^∘ E) when assessing trends. For SW↓ and LW↓ only the station in Bergen (UiB) meets the criteria. For consistency the humidity observations from Bergen-Florida are inspected for trends as well.

6 Comparison between existing long-term estimates and the new hybrid approach, HySN

Daily estimates of near-surface humidity and SW↓ and LW↓ from MERRA2, ERA-Interim, PGMFDv2, WFDEI, VFDv1, VFDv2, and HySN (see Table 1) from 1982 to 1999 are first inspected for vertical gradients in order to adjust the model estimates to station altitude in the following. Further, the estimates' quality on a multi-annual timescale is assessed by considering their mean and absolute deviations from station measurements. The model estimate's distribution is assessed by comparing their daily cumulative distribution to that of the station observations. The estimate's similarity to the observations on a daily timescale is considered by inspecting the anomaly correlation coefficient, i.e., their daily correlation with the measurements after the seasonal cycle has been subtracted, and by considering if their differences to station measurements show sensitivity to the local daily air mass type, which has been classified for seven Norwegian stations by Bower et al. (2007). For humidity and SW↓, the geographical gradients in the models are compared with those in the observations using multiple linear regression. A separate subsection presents modeled and observed trends in each calendar month from 1985 to 1999. Humidity trends are compared after grouping the observations and model estimates into five regions, while SW↓ and LW↓ are computed for Bergen, the only location where long-term measurements of SW↓ and LW↓ are available within the time period with little missing data.

6.1 Humidity

The model estimates of near-surface humidity are compared against humidity observations from 84 stations. The observation stations used in the validation of humidity are located on average 200 m below the coarse-scale grid cell elevation, i.e., in the fjords and on the coast rather than in the surrounding terrain.

6.1.1 Vertical gradients in humidity

Multiple linear regression of seasonal mean humidity at the location of the humidity stations shows that altitude is a significant predictor of humidity in the observations and all models (see Sect. S2 of the Supplement). The vertical gradient in the observations is close to the moist adiabatic lapse rate but varies considerably with season and distance to the coast (C). On average, dew point temperature decreases by 5.2 ^∘C km⁻¹ increase in altitude in summer and freeze point temperature decreases by 4.4 ^∘C km⁻¹in winter. Regression based on vapor pressure is found to give smaller relative errors than regression based on dew point temperature. This is because dew point temperature has a higher sensitivity to temperature at low temperatures. The observed vertical gradient in vapor pressure is −0.39 hPa per 100 m in summer and −0.24 hPa per 100 m in winter, and the vertical gradient is weakened by 0.11 hPa per 100 m every 100 km away from the coast.

The vertical gradients in the humidity data differ depending on the data source; e.g., the estimates from PGMFDv2 and WATCH show a weaker decrease with altitude than the observations and other models. For each model the vertical gradients are computed as the mean of the seasonal coefficients found in the regression analyses of the model in question and the observations, linearly interpolated from a seasonal to a daily frequency. The altitudinal adjustment results in a mean difference in humidity, expressed as dew point temperature (T_d) of about 1 ^∘C for the coarse-scale models and about 0.06 ^∘C for the estimates with a 1 km×1 km resolution. The largest adjustment is an increase in MERRA2's T_d of 7.3 ^∘C at Tafjord station, where the model's orography is 1154 m above the station altitude.

6.1.2 Differences of humidity estimates to station observations

The seasonal cycle of the observations and models (adjusted to station altitude) is shown in Fig. 3a. The largest deviations are seen in the period of highest humidity, i.e., during summer. The signs of the average deviations are consistent throughout the year. PGMFDv2 (denoted with P), MERRA2 (M), and to some degree ERA-Interim (E) and WFDEI (W) show larger estimates than the observations, whereas VFDv1 (V1) and VFDv2 (V2) generally show lower estimates. The HySN (H) estimates follow the mean monthly values of the observations closely. This is also evident in Table 3, where summary of statistics for the humidity estimates presented, and HySN shows a mean station error in T_d of just 0.1 ^∘C. An aggregated mean similar to the observations does not ensure small deviations from the measurements at individual stations.The VFD estimates have the second smallest deviation in aggregated mean (Δ); however, when considering the average absolute deviation ($\left|\mathrm{\Delta }\right|$) HySN, WFDEI, and ERA-Interim perform better than VFD.

Figure 3 The seasonal cycles of monthly 2 m vapor pressure (a), incident shortwave radiation (b), and incident longwave radiation (c), averaged over the location observations are available. The month of the year is denoted on the horizontal axis. The observations (O) are plotted with a thick, continuous black line. MERRA2 (M, solid) and ERA-Interim (E, dashed) are plotted in red, PGMFDv2 (P, solid) and WFDEI (W, dashed) in blue, VFDv1 (V1, solid) and VFDv2 (V2, dashed) in orange, and HySN (H) with a dashed lilac line.

Table 3 Summary of metrics showing the humidity estimates' similarity to station observations. Differences (Δ) are given in dew point temperature in degrees Celsius. Δ is the mean station difference, $\left|\mathrm{\Delta }\right|$ is the mean absolute station difference, $\left|{\mathit{\delta }}_{max}\right|$ is the largest absolute difference at any station, while $|{\mathit{\delta }}^{\mathrm{s}}{|}_{max}$ is the largest seasonal difference at any station. ACC is the anomaly (de-seasonalized) daily correlation coefficient, while K-S indicates the number of stations where the daily mean cumulative distribution passes the Kolmogrov–Smirnov test of similarity (p>0.001). The best scores are shown in bold.

Differences between the model estimates and the station measurements of humidity, expressed as dew point temperature, are depicted for each station, sorted from south (upper y axis) to north (lower y axis) in Fig. 4. The mean absolute difference ($\left|\mathrm{\Delta }\right|$ or MAE) varies from 0.7 ^∘C for the HySN estimates to 1.8 ^∘C for the PGMFDv2 estimates. The largest deviation occurs at an inland station, Fagernes, where PGMFDv2 T_d estimates a 5.4 ^∘C higher T_d than observed. The figure suggests a latitudinally dependent bias for certain models, and this is further explored in the following subsection by evaluating the models' geographical gradients.

Figure 4 For each station, sorted from south to north (y axis), and each model (x axis) the differences between the modeled and observed station mean dew point temperature (a), incident shortwave radiation (b), and incident longwave radiation (c) are shown.

The humidity estimates are evaluated on a daily basis by de-seasonalizing the time series (subtracting the observed day-of-year mean). Figure 5 shows, for each model, the de-seasonalized time series of humidity, the mean temporal correlation coefficient (now equivalent to the anomaly correlation coefficient, ACC), the mean normalized root-mean-square error, and the mean standard deviation in a Taylor plot. Figure 5a visualizes the mean station metrics for the de-seasonalized time series of humidity. The estimates from ERA-Interim and post-processed ERA-Interim (HySN and WFDEI) are closest to the observations and show similar results. MERRA2 also shows a high ACC. PGMFDv2, which is based on an older reanalysis with lower spatial resolution, and the VIC type estimates show slightly poorer results, with an ACC ranging between 0.5 and 0.7 (see also Table 3).

Figure 5 Taylor plots depicting the standard deviation ratio and correlation coefficients (ACCs) for the de-seasonalized time series of vapor pressure (a), incident shortwave radiation (b), and incident longwave radiation (c).

6.1.3 Evaluation of geographical humidity gradients

Multiple linear regression models are fitted to seasonal mean humidity with the four seasons as categorical predictors, where fall (autumn) is the baseline season in the model. The geographical predictors considered are altitude (z given in kilometers), latitude (above 57^∘ N), and distance to the coast (C, given in per 100 km increments). Further, interactions between altitude and season and between altitude and continentality are included. Each model is paired with the observational data in a common regression model where the data source is included as a categorical predictor.

Figure 6 displays the regression coefficients for the observations and the coefficients for the models if they are significantly different (p<0.01) from those of the observations. Higher significance is marked with a darker color. The HySN estimates have similar coefficients to the observations. The regression shows (Fig. 6) that MERRA2 and PGMFDv2 have significantly higher intercepts (higher fall mean values at the coast of southern Norway) than the observations. MERRA2 further shows a stronger latitudinal gradient and a much weaker decrease in humidity with distance from the coast than the observations. In addition to having a higher intercept than the observations, PGMFDv2 shows a more pronounced seasonal dependency, a weaker continental gradient, and a 50 % stronger latitudinal gradient than the observations. VFDv1 and VFDv2 show a more than 60 % more pronounced decrease in humidity with continentality than the observations. VFDv2 also shows a weaker increase in humidity in summer than the observations.

Figure 6 Seasonal and geographical dependencies of seasonal humidity (a) and daily clearness index (b) are depicted. The row names are the names of the coefficients, including the intercept (I), of the multiple linear regression model. The regression coefficients of the observational data are shown in the leftmost column of each plot, while the coefficients found for the model estimates are only shown if they are significantly different from those of the observations (using a limit of p<0.01 for humidity and p<0.05 for CI). Lower p values are indicated with darker colors using a logarithmic color scale (log ₁₀(p)).

6.1.4 Air mass type sensitivity of humidity deviations

In Fig. 7a the daily deviations of the humidity estimates are grouped according to the location's daily air mass type classification (SSC type; see Sect. 5). The classification is available for Sola, Fornebu, Flesland, Bodø, Tromsø, and Slettnes stations and the humidity observations considered are from Saerheim, Aas, Bergen, Bodø, Tromsø, and Kirkenes stations. All the estimates are too humid in dry weather types. The PGMFDv2 and VFDv2 estimates show considerable overestimations of humidity in dry weather types and underestimations in moist weather types. The lack of range is consistent with the lower normalized standard deviation seen in the Taylor plot (Fig. 5a). The ERA-Interim, WFDEI, and HySN differences also show a slight sensitivity to air mass type but much less than the VFDv1, VFDv2, MERRA2, and PGMFDv2.

Figure 7 Differences between model estimates and observed values binned according to the daily air mass type, classified for nearby stations within Norway (Bower et al., 2007). The air mass types are dry polar (DP), dry moderate (DM), dry tropical (DT), moist polar (MP), moist moderate (MM), moist tropical (MT), and transitional (T). Panel (a) shows the mean difference in vapor pressure for Saerheim, Aas, Bergen, Bodø, Tromsø, and Maze stations. Panel (b) shows differences between model estimates and observed incident shortwave radiation at Saerheim (Sola airport), Aas (Fornebu airport), Bergen (Flesland airport), and Trondheim (Ørland airport). Panel (c) depicts differences in incident longwave radiation binned according to air mass type in Bergen (Flesland airport) and Trondheim (Ørland airport).

6.2 Incident global shortwave radiation (SW↓)

SW↓ observations from 10 sites on the Norwegian mainland are considered. At most locations the coarse-scale models' corresponding grid cells have an altitude 300–400 m above station altitudes.

6.2.1 Vertical gradients in clearness index (CI)

Multiple linear regression was used to provide a vertical gradient in SW↓, expressed as clearness index (CI, i.e., the fraction of SW↓ of the extraterrestrial incoming radiation, SW↓_E), in order to adjust the estimates to the stations' altitudes. Multiple linear regressions, including both continentality and altitude as predictors resulted in altitudinal coefficients varying in both magnitude and sign for the different models and observations (not shown). This was likely because the correlation between altitude and continentality varies between 0.56 and 0.86 depending on the data source. Excluding continentality from the predictors provided vertical CI gradients with a consistent sign. The observations show vertical CI gradients of 0.020∕100 m in winter, 0.013∕100 m in spring, 0.005∕100 m in summer and 0.003∕100 m in fall (see Fig. 6b). The observed SW↓ thus increases, on average, with altitude in all seasons.

The effect of adjusting the model estimates to station altitude is an average reduction in SW↓ of 0.7–1.5 W m⁻² for the coarse-scale models (MERRA2, ERA-Interim, PGMFDv2, and WFDEI) and a reduction of merely 0.1–0.3 W m⁻² for the models with a 1 km×1 km grid (VFDv1, VFDv2, HySN). The largest adjustment is a mean reduction of the PGMFDv2 SW↓ estimate of 4 W m⁻² at a station in southeastern Norway (Gjengedal).

6.2.2 Differences of SW↓ estimates to station observations

The mean monthly model estimates of SW↓ averaged over all 10 stations, after adjustment to station altitude, are visualized in Fig. 3b. In winter the deviations are small, but in spring and summer all models except VFDv1 overestimate SW↓. The MERRA2 SW↓ is on average 35 W m⁻² higher than the observations in July, and the VFDv2 SW↓ is 29 W m⁻² higher than the observations in both June and July. ERA-Interim, WFDEI, and HySN show the largest overestimations in May, a month when solar radiation is high and snow cover is variable.

Figure 4b depicts the mean difference between the model estimate and the observations of SW↓ at individual stations. At half of the stations the mean difference between the WFDEI and HySN and the observations is lower than the measurement uncertainty of newer pyranometers in optimal conditions (Sect. 4). The figure further shows that most models consistently overestimate SW↓. This is not true for VFDv1. While VFDv1 has the second lowest mean monthly deviation (Fig. 3), its mean absolute difference is larger, on average 10 W m⁻². This is also evident in Table 4 where summary statistics for the SW↓ estimates are presented.

Table 4 As in Table 3 but with metrics listed for SW↓. Except for ACC and K-S, which are dimensionless, the units are given in W m⁻².

Table 4 shows that at individual stations seasonal deviations in model estimates from station observations are as large as −62 W m⁻². The large underestimation is found in VFDv1 and not VFDv2 at a coastal station in northern Norway, Bodø. Also listed in the table is the percentage of stations where the daily model estimate, adjusted to station altitude, passes the Kolmogorov–Smirnov test of similarity of their cumulative distribution with the observations, which is zero cases for PGMFDv2 and 70 % for HySN.

The similarity of the model estimates to the observations at a daily frequency is visualized in a Taylor plot (Fig. 5b). As also seen for the humidity estimates, PGMFDv2, VFDv1, and VFDv2 have lower ACCs (31 %–48 %) than the estimates based on newer reanalysis data (73 %–78 %). PGMFDv2 in particular shows a variance at a daily frequency that is considerably smaller than the observations.

Figure 8 Mean daily clearness index (CI) during winter (a) and summer (b) is depicted in violin plots, where the kernel density distribution of the observations and each model is shown, mirrored across the y axis. The observed median value is drawn with a blue solid line.

6.2.3 Evaluation of geographical gradients in clearness index

Similar to what has been done for humidity, the observations and the corresponding vertically adjusted model estimates of daily CI are compared using multiple linear regression. The seasons, latitude, and altitude are used as predictors, including interaction between season and latitude and between season and altitude. Figure 6b shows the regression coefficients of the observational data in the leftmost column. The coefficients found for the model estimates are only displayed if they are significantly different (p<0.05) from those of the observations. Larger differences (log ₁₀(p)) are marked with a darker color. The estimated intercept of PGMFDv2 stands out in the plot. It is 40 % higher than the estimated intercept of the observations. The second most evident difference is the latitudinal gradient in both VFDv1 and VFDv2, which is several times stronger than observed. PGMFDv2 also shows a stronger latitudinal gradient than observed. Other notable differences are the estimated summer CI values of MERRA2, which are considerably higher than those seen in the observations.

6.2.4 Air mass type sensitivity of SW↓ deviations

Figure 7b shows the differences between the models' SW↓ estimates and observations at Aas, Saerheim, Bergen-GFI, Bodø, and Tromsø stations, grouped according to the weather type at nearby weather stations (Fornebu, Sola, Flesland, Bodø, and Tromsø, respectively). All models except VFDv1 show positive SW↓ deviations during weather types classified as moist, which occur most frequently. In the less prevalent dry weather types all models except MERRA2 show slightly lower estimates than observed. The largest underestimations are seen for the VFDv1 estimates. Further, the deviations of PGMFDv2, VFDv1, and VFDv2 show a stronger dependency on weather type than MERRA2, ERA-Interim, WFDEI, and HySN. The MERRA2 SW↓ estimates are, however, overestimated during all weather types. The WFDEI SW↓ estimates show considerable deviations from the ERA-Interim estimates, with larger underestimations found during both dry polar and dry tropical weather types, and considerably lower overestimations found on days classified with a moist tropical weather type. Grouping the clearness index into either dry or moist and transitional weather types shows that the observed CI decreases on average by 0.22 in moist and transitional types. A similar decrease is seen in MERRA2, ERA-Interim, WFDEI, and HySN but not in VFDv1 and VFDv2 (−0.12) or PGMFDv2 (−0.05). The summer and winter distributions of clearness index at the stations considered are depicted in Fig. 8. It is evident that PGMFDv2 spans a much smaller range of transmissivity than observed in both summer and winter and that the VIC type estimates have a bias towards low estimates and show less variability than observed in winter.

Since both WFDEI and HySN are based on ERA-Interim, and ERA-Interim shows overestimations of SW↓ in summer, where observations were available, the differences in estimated SW↓ in ERA-Interim and the observations were inspected for dependencies on differences in modeled and observed cloud cover, near-surface humidity, and snow cover using regression. The results varied with both season and location but for the aggregated data significant dependence on differences in observed and modeled 2 m humidity and snow cover was found (with higher snow cover in the model associated with higher SW↓ estimates in the model), and in the warm season larger overestimations were seen in ERA-Interim when the model produced high clouds.

6.3 Incident longwave radiation (LW↓)

Only two stations have LW↓ observations available during the validation period, Bergen-GFI (western Norway) and Trondheim-Voll (central Norway); both are located near the coast (see Fig. 1 and the Supplement). More than 17 years of daily measurements are available from the Bergen-GFI station, while at Trondheim-Voll only about 2 years of observations are available.

6.3.1 Vertical gradients in LW↓

Since only two stations have longwave observations, no altitudinal gradient can be inferred from the observations; instead altitudinal gradients are taken from ERA-Interim. The previous comparison of altitudinal gradients within the observations and models has shown that ERA-Interim has similar altitudinal gradients to the observations. Further, ERA-Interim has not previously been vertically adjusted. The vertical gradients found ranged from −4.0 W m⁻² per 100 m in December to −0.6 W m⁻² per 100 m in June, were weakened by 0.20 W m⁻² per 100 m for every 10 km away from the coast, and strengthened by 0.23 W m⁻² per 100 m for every latitude north of 57^∘ N. On average the vertical gradients were around −4.5 W m⁻² per 100 m in winter and −1.8 W m⁻² per 100 m in summer. The gradients were temporally interpolated to day-of-year values and applied to the models in order to adjust the estimates to the altitude of the two observational stations. The largest change in the estimates due to altitudinal adjustment is an average increase of 8.6 W m⁻² for MERRA2 when adjusted to the 278 m lower altitude of Bergen station compared to MERRA2's grid cell altitude.

6.3.2 Differences of the LW↓ estimates to station observations

Figure 3c depicts the mean monthly LW↓ at the two stations. Summary statistics for the LW↓ estimates after adjustment to station altitude are also presented in Table 5. At the two stations all models except WFDEI estimate lower values than observed in all months. WFDEI, however, estimates on average 10–15 W m⁻² more LW↓ than is observed from October through January. The largest absolute differences are found in MERRA2, where LW↓ is underestimated at 8 W m⁻² in winter and at 25 W m⁻² midsummer. The remaining models estimate 11 to 17 W m⁻² less LW↓ than observed in summer, and show smaller underestimations in winter.

Table 5 As in Table 3 but with metrics listed for LW↓. Except for ACC and KS, which are dimensionless, the units are in W m⁻².

The skill of the model estimates in capturing the day-to-day variability in LW↓ is visualized in Fig. 5c, indicating the correlation and normalized standard deviation of the de-seasonalized time series. The estimates based on newer reanalysis data, MERRA2, ERA-Interim, WFDEI, and HySN have anomaly correlation coefficients around 0.8, while PGMFDv2, VFDv1, and VFDv2 have lower ACCs (46 %–60 %).

6.3.3 Air mass type sensitivity of LW↓ deviations

Figure 7c depicts the differences between the model's LW↓ estimates and the observed values at Bergen and Trondheim stations, grouped according to the daily weather type (classified at Flesland and Ørlandet stations). On days classified with moist or transitional weather types, all models except WFDEI underestimate LW↓. PGMFDv2, VFDv1, and VFDv2 clearly have weather-type-dependent deviations from the observations, with underestimations in moist weather types and smaller underestimations or even overestimations compared to the observations in dry weather types. The MERRA2, ERA-Interim, WFDEI, and HySN estimates largely show similar differences to the observed values in all weather types. An additional comparison of the ERA-Interim estimates at Bergen, where cloud observations are available, showed no difference in average deviation in the estimates of incident longwave radiation on common clear-sky days compared to the remaining days (not shown). The lower estimates of incident longwave radiation in ERA-Interim are thus not likely to be primarily related to differences in cloud properties.

6.4 Modeled and observed trends in near-surface humidity, SW↓, and LW↓ from 1985 to 1999

January 1985 is considered to be the start of the brightening period in Europe after a period of SW↓ dimming due to aerosol emissions (see, e.g., Wild et al., 2005). In the following, where available, observations and co-located model estimates are inspected for trends in near-surface humidity, SW↓, and LW↓ from 1985 to 1999.

After screening the observational time series 59 humidity stations are grouped into five geographical regions (southwest (S_W), southeast (S_E), central (C), northwest (N_W), and northeast (N_E); see Sect. 5). Table 6 lists the results of the trend tests for each calendar month and region, listing only the Sen slope if the Mann–Kendall trend test is significant at a 5 % or 10 % level, with the latter denoted in italics. In all regions except the northeastern part of Norway significant increases in observed T_d occurred in September from 1985 to 1999. The observations further show an increase in T_d in April in southern Norway, an increase July in the northeastern part of Norway, and a decrease in May T_d in central Norway. For southern and central Norway all of the models capture the increase in September T_d. The models reproduce the observed trends in spring and in northern Norway to a lesser degree.

Table 6 Linear, decadal trends in monthly mean T_d (^∘C) between 1985 and 1999, significant at a 5 % or 10 % (denoted in italics) level in the observations (O) and the model estimates, listed with the month and region denoted on the left (month, region).

Bergen is the only location in Norway where long-term records of incident shortwave or longwave radiation are available with little missing data within the time range. Between January 1985 and December 1999, observations from the University in Bergen show trends in annual SW↓ of 1.7 W m⁻² per decade (p<0.1), while LW↓ decreases with −8.4 W m⁻² per decade (p<0.001). At the nearly co-located measurement station of the Norwegian Meteorological Institute, at Bergen-Florida, annual dew point temperature has a trend of 1.2 ^∘C per decade (p<0.0001).

In individual calendar months larger trends were found. Observed August mean SW↓ increased with 51 W m⁻² per decade (Table 7). The observations also show modest, but significant, increases in October and December. Apart from WFDEI and VFDv1, which show no significant SW↓ trends, the models reproduce a significant increase in SW↓ in August, and no models show trends of the opposite sign during the time period considered.

Monthly mean LW↓ in Bergen shows a significant decrease during several months of the year: the largest is found in May, with −21 W m⁻²per decade (Table 8). In August, October, and December, months when concurrent increases in shortwave radiation were found (Table 7), the observations show significant decreases in incident longwave radiation. None of the models show equally strong negative trends in monthly mean LW↓ in their corresponding grid cells. MERRA2 and ERA-Interim show no significant trends during the time period, whereas PGMFDv2, WFDEI, and HySN exhibit one single month with a negative trend. Meanwhile, VFDv1 and VFDv2 show weakly positive trends in September. The increasing trend in incident longwave radiation in September in the VIC type estimates may be related to the concurrent increase in humidity. All the models also show an increasing trend in humidity in the grid cell covering Bergen in September; however, the models generally show weaker trends than observed at Bergen-Florida station (see Table 9).

Table 7 Linear, decadal changes in monthly mean SW↓ (W m⁻²) at Bergen between January 1985 and December 1999, significant at a 5 % or 10 % (denoted in italics) level.

Table 8 Linear, decadal changes in monthly mean LW↓ (W m⁻²) at Bergen between January 1985 and December 1999 significant at a 5 % or 10 % (denoted in italics) level.

Table 9 Linear, decadal changes in monthly mean T_d (^∘C) for Bergen-Florida station and for the co-located grid cells of the model estimates between January 1985 and December 1999, significant at a 5 % or 10 % (italics) level.

7 Discussion

Historical estimates of humidity and incident shortwave and longwave radiation have been compared to station observations from mainland Norway from 1982 through 1999. A total of 84 stations provide vapor pressure (VP) observations, 9 stations provided SW↓ observations, while only 2 stations provided LW↓ observations. The estimates evaluated are from two reanalysis data sets, MERRA2 and ERA-Interim; three data sets composed of reanalysis data blended with gridded observational data, PGMFDv2, WFDEI, and HySN; and two versions of the VIC type forcing data, estimates based on gridded observational data combined with empirical algorithms.

Differences between the estimates and observations are not necessarily due to errors in the estimates, as a vertical adjustment to station altitude is not a sufficient reason to require that the model grid cell estimates should equate to the observed values. The numerical model estimates may still differ from the observations for valid reasons, such as differences in snow cover, differences in land cover type (the observations are from sensors usually located over grass or, in some cases, on top of buildings), and the averaging out of sub-grid variability in the models (see, e.g., Göber et al., 2008). The uncertainty in the observations may also contribute to the differences. However, large differences may suggest biases in the estimates.

7.1 Vertical gradients

Significant vertical gradients were found for humidity, incident shortwave, and incident longwave radiation, justifying an altitudinal adjustment to station altitude before comparison of the model estimates with the station observations (ℋ_a). The altitudinal vapor pressure gradients found here were on average −0.25 hPa per 100 m in winter and −0.34 hPa per 100 m in summer. The summer gradient is similar to what Marty (2000) found in the Alps in summer, −0.35 hPa per 100 m; however, the winter gradient is considerably higher than found in the Alps (−0.14 hPa per 100 m). The impact of adjustment to station height was small for the estimates with a finer spatial scale, only a 0.06 ^∘C change in T_d on average, while for the coarser-scale estimates, MERRA2, ERA-Interim, PGMFDv2, and WFDEI, the impact of the vertical adjustment was considerably larger, resulting in an average 1 ^∘C increase in T_d. The WFDEI and PGMFDv2 showed weaker vertical humidity gradients than observed. This may be a result of the interpolation techniques employed in the CRU T₂ data set used to bias-correct and downscale both the ERA-Interim and NCEP-NCAR Reanalysis, or due to the use of a constant temperature lapse rate of 6.5 ^∘C when interpolating the temperature of the two reanalyses to the CRU orography. Notably, the vertical gradients in near-surface humidity in MERRA2, a reanalysis where surface observations are not assimilated (see Table 1), are similar to the vertical gradients found in the observations and those found in ERA-Interim.

Observed SW↓ in the form of clearness index (CI; see Sect. 5) showed the highest altitudinal gradient in winter, a slightly lower gradient in spring, and rather low gradient in summer and fall. The vertical gradients found are larger than the gradient of 0.00295 per 100 m used in the implementation of the Bristow and Campbell (1984) algorithm in historical versions of the VIC preprocessor (MTCLIM versions before 4.2, before the Thornton and Running, 1999, algorithm was implemented). Though the CI gradient is stronger in winter, the considerably smaller amount of SW↓ received leads to a weaker gradient in SW↓. The CI gradients translates to SW↓ gradients of about 0.3 W m⁻² per 100 m in fall and winter, 1.6 W m⁻² per 100 m in spring, and 1.2 W m⁻² per 100 m in summer. Marty (2000) found all-sky gradients in SW↓ in the Alps of 1.1 W m⁻² per 100 m in winter and 0.7 W m⁻² per 100 m in summer. The differences between the gradients found here and those given in Marty (2000) may likely be explained by the differences in the received extraterrestrial radiation and differences in cloud and snow cover climatology. The models largely showed similar vertical CI gradients to the observations. The exceptions were PGFMDv2 and VFDv1; PGMFDv2 showed significantly (p<0.01) weaker vertical gradients with a weaker seasonality, and VFDv1 produced a stronger vertical CI gradient in summer than in winter. The adjustment of the coarser-scale estimates resulted, on average, in a 5 times larger change in SW↓ for the coarse-scale estimates than the finer-scale estimates. The regression-based vertical adjustment produced similar SW↓ estimates for ERA-Interim and the HySN estimates.

Since two LW↓ stations located more than 400 km apart could not provide an observation-based vertical gradient in LW↓, ERA-Interim was consulted instead. The gradients were, on average, −4.5 W m⁻² per 100 m in winter and −1.8 Wm⁻² per 100 m in summer. Marty (2000) found vertical gradients in incident longwave radiation of −2.8 Wm⁻² per 100 m in winter and −3.1 Wm⁻² per 100 m in summer for the Alps and of −4.1 W m⁻² per 100 m in winter and −2.6 W m⁻² per 100 m in summer, when considering a subset of observation stations in Switzerland. The different vertical gradients found may be explained by differences in temperature and humidity gradients, different climatological distributions of clouds, and the difference in initial temperature, as LW↓ is a function of temperature to the power of 4. The regression-based vertical adjustment of ERA-Interim LW↓ estimates resulted in a larger correction of LW↓ than the clear-sky adjustment implemented in HySN, alluding to the fact that the clear-sky altitudinal adjustment implemented in similar data products might be too low, especially for locations with a maritime climate, like Bergen.

7.2 Differences to station observations

7.2.1 Humidity estimates

The empirically based model estimates, VFDv1 and VFDv2, show, on average, slightly lower estimates of humidity than observed. Both VFD type estimates are found to show a 50 % stronger decrease in humidity with continentality than the observations (see Sect. 6.1.3). The modified version of the Magnus type formula, based on Kimball et al. (1997), used in MTCLIM to generate the VFD humidity estimates is likely not appropriate for Norway. Previous studies, e.g., in the development of gridded climate variables by New et al. (1999) and in the application of the MTCLIM model over complex terrain in Australia (Thornton et al., 2000) and in the western US (Pierce et al., 2013), found that the Kimball et al. (1997) method did not result in overall improved humidity estimates. Indeed, in Kimball et al. (1997) the method is found to give improved estimates of humidity in locations where the ratio of potential evaporation to annual precipitation is larger than 2.5. In most regions of Norway this ratio is well below unity. The more conventional method of using daily minimum temperature as a proxy for dew point temperature will likely give relatively small overestimations of humidity compared to the underestimations resulting from using the Kimball et al. (1997) method.

The reanalysis-based estimates all overestimate humidity, and the overestimations are generally higher in weather types classified as dry according to the methodology of Bower et al. (2007). MERRA2 and PGFMDv2 particularly overestimate humidity in dry conditions. The same two models also show a significantly stronger decrease in humidity with latitude than observed. MERRA2 also shows a weaker decrease in humidity with continentality. The weaker decrease in humidity with continentality seen in MERRA2 may perhaps be partly explained by the model's coarse resolution and land mask (see Fig. 1), and MERRA2's exaggerated latitudinal gradient in humidity in Norway may perhaps be associated with MERRA2's larger latitudinal gradient in SW↓.

The humidity estimates from HySN match the observations best, considering all metrics except from the anomaly correlation coefficient (ACC). The ACC of ERA-Interim estimates is marginally higher (0.02) than in HySN. This is likely due to the capping of relative humidity at 100 % when applying the SeNorge2 temperature in the development of HySN. Combining the methods outlined in Cosgrove (2003), which for humidity relies on the assumption of constant relative humidity with altitude, a high-quality reanalysis data set (ERA-Interim), and a high-resolution, national, observation-based temperature data set is found to provide high-quality daily estimates of humidity in the current study region. The coarser, reanalysis-based data sets generally show higher ACCs than the VFD estimates. Numerical weather prediction (NWP) models are skilled at capturing synoptic events, i.e., weather or climatological patterns on a spatial order of 1000 km, and a temporal order of days or weeks, such as cold air outbreaks and the changing sources of air masses during the passage of warm and cold fronts. Though the NWPs may have systematic biases and a much lower spatial resolution than empirically based estimates, it is not surprising that they are useful in representing daily weather variability.

7.2.2 Incident shortwave radiation

Shortwave incident radiation is, on average, overestimated for all model estimates except VFDv1. HySN, ERA-Interim, and WFDEI vary in obtaining the highest ranking depending on the metric considered. For instance, WFDEI shows a slightly lower average deviation from the observations than ERA-Interim and HySN. On the other hand, WFDEI shows larger underestimations in dry weather types than ERA-Interim and HySN (Fig. 7). Overall, the three models provide vertically adjusted estimates of incident SW↓ close to the observations, with average deviations from station measurements below 4 W m⁻² and ACCs above 0.76.

The average difference between the ERA-Interim estimates and the observations is smaller than in Urraca et al. (2018), where an average overestimation of 12 W m⁻² was found when comparing ERA-Interim SW↓ estimates to station measurements in Europe between 2010 and 2014. The smaller difference seen in the current study may in part be explained by the relatively small amount of solar radiation reaching Norway, the different time periods considered, and the vertical adjustment included in the current study. Urraca et al. (2018) also found that MERRA2 shows poorer results than ERA-Interim, with average overestimations of 18 W m⁻². This is consistent with our findings, where MERRA2 has the highest mean deviation from the observations of any of the considered estimates. Overestimations of incident shortwave radiation over land are not only an issue of reanalysis data sets covering Europe but have been a long-standing issue in global (Wild et al., 2015) and regional climate models (Katragkou et al., 2015; Jerez et al., 2015).

Two versions of VIC type forcing data are evaluated in the current study. The two versions differ in their input data and in the version of VIC preprocessor used. The oldest version of the VFD data sets is partly based on a 11 km national reanalysis (NORA10) to provide maximum and minimum temperature. The older version showed large underestimations of incident shortwave radiation at several stations, particularly near the coast in northern Norway (Fig. 4). These findings are in line with Bohn et al. (2013), where the MTCLIM algorithms were found to underestimate incident SW↓ radiation by 26 %, on average, at coastal sites. The MTCLIM algorithms implemented in VFD rely in part on the diurnal temperature range to estimate cloud cover, using a low range as an indication of cloud cover. Near the coast, the diurnal temperature range may be low due to the moderating influence of the nearby ocean, due to its high heat capacity. The more recently compiled version of VFD data, VFDv2, which is based on a newly developed, high-resolution gridded data set of T_min and T_max, does not show similar underestimations near the coast of northern Norway. The different estimates produced indicate that great care must be taken to make sure the VIC style forcing data have consistent input data and algorithm versions if the data are used in, for example, climate change impact studies.

The newer and higher-resolution input data used in VFDv2 did not result in a lower mean absolute station deviation, as its SW↓ estimates were consistently overestimated. Both VFD versions show a much stronger latitudinal gradient than observed and a too strong altitudinal gradient in summer. The latter finding is in line with Mizukami et al. (2014), where VFD type estimates for the Colorado River basin showed increasing overestimations of SW↓ with increasing altitude. The exaggerated latitudinal gradient in SW↓ may be connected to the use of the diurnal temperature range in the algorithm. Bohn et al. (2013) found that the relationship between cloud cover and the diurnal temperature range breaks down for ranges below 5 ^∘C. Further, New et al. (1999) states that the relationship between diurnal temperature range and cloud cover is weak at around 60^∘ N in winter and, further, becomes positive in the Arctic.

Binning the estimates according to air mass type shows that the PGFv2 and VIC type estimates show less sensitivity to the prevailing weather type than the observations. The observations and the remaining model estimates show a decrease in clearness index of about 0.22 on days classified as moist or transitional weather types rather than dry, while the VIC type estimates and PGFv2 show reductions of 0.12 and 0.05, respectively. On average, the VIC type estimates and PGFv2 underestimate incident radiation in dry weather types (see Fig. 7). The similarity between the PGMFDv2 and VFD estimates of SW↓ may be explained by the fact that the PGMFD SW↓ is bias-corrected based on gridded cloud cover from CRU using the Thornton and Running (1999) relationship between SW↓ and cloud cover, which is also used in VFD. Further, the gridded CRU cloud cover data set is a secondary or derived observational data set, which is, similar to VFD, in part based on regression using diurnal temperature range as a predictor (New et al., 1999). The lower sensitivity to air mass type found in PGMFDv2 and the VIC type forcing data might contribute to the lower ACC found for these estimates.

7.2.3 Incident longwave radiation

The evaluation of incident longwave radiation is compromised by the lack of observational data. Only two sites observe incident longwave radiation in the considered time period. The difference between the annual mean of the model estimates and observations are for the two stations considered larger for incident longwave radiation than for incident shortwave radiation. The annual deviations ranges from −16 to +7 W m⁻². Underestimations of monthly means are found throughout the year for all models except WFDEI. The deviations from the station observations for WFDEI, MERRA2, ERA-Interim, and HySN are largely similar in all weather types; i.e., the underestimations are also found on days classified with dry weather types. An additional evaluation of the ERA-Interim LW estimates for Bergen, where cloud observations are available, showed that the deviations from station observations were similar on days where clouds were present in the observations, the model, and the remaining days.

While overestimation of incident shortwave radiation has been a long-standing issue in many climate models and reanalyses, incident longwave radiation is typically underestimated (Katragkou et al., 2015; Li et al., 2016; Zib et al., 2012; Wild et al., 2017). The causes of the underestimation are, however, debated. Li et al. (2016) points to an improper representation of interactions between radiation and suspended frozen water particles in the atmosphere (solid hydrometeors) as a culprit, while Zib et al. (2012) speculate that errors in simulated aerosols, water vapor content, and cloud properties (rather than cloud amounts) are the cause. Local issues such as longwave emissions from nearby terrain may also contribute to the deviations (Rontu et al., 2016). Lastly, observational uncertainty confounds the picture further, particularly given that the two sensors were unshaded.

The anomaly correlation coefficients are, as also seen for humidity and incident shortwave radiation, considerably lower for PGFMDv2 and the VFD estimates than the estimates based on newer reanalysis data (0.46–0.60 vs. 0.79–0.82). This may be caused by the representation of clouds in the models. As discussed for the PGFMDv2 and the VFD SW↓ estimates, the use of diurnal temperature range as a proxy for cloud cover may not be suitable for the current maritime, high-latitude study region.

7.3 Trends

The analysis of observed humidity trends from 1985 to 1999 showed significant increases in April in southern Norway, a decrease in May in central regions of Norway, significant increases in July in the northeastern part of Norway, and increases in all regions except the northeastern part of Norway in September. All the data sets, both the reanalysis-based estimates and the more empirically based VFD estimates, capture the increase in humidity seen in September. The significant increases found in humidity when averaging over the stations in southeastern and southwestern Norway in April are not seen in any of the models. The VFD estimates do, however, capture some of the increases in humidity that were seen in the measurements from Bergen-Florida (Table 9). A recent study by Nilsen et al. (2017) found that changes in large-scale weather patterns can, in part, explain the significant increases in 2 m temperature between 1981 and 2010 seen in Scandinavia in September but not most of the increases seen in April. Another inquiry by Rizzi et al. (2017) found that the increasing temperature trends seen in May in many parts of Norway showed a strong correlation with a concurrent decrease in snow cover. The decline in snow cover in May found in Rizzi et al. (2017) was particularly strong in low-lying areas. If the changes in temperature and humidity are connected to local changes in snow cover, it is possible that the coarser-scale reanalysis data, which often have a mean grid cell altitude above the measurement station elevation, do not capture the measured changes.

Surface incident radiation was inspected for trends from 1985 to 1999 at the one station where measurements are available in the time period with little missing data: Bergen. A hardly significant (p<0.1) annual trend in SW↓ was found in the observations, 1.7 W m⁻² per decade. However, in individual calendar months larger trends were found. The largest trend, 51 W m⁻² per decade, was found in the observations in May. The observed August trend was reproduced fairly well in ERA-Interim, PGMFDv2, and HySN, and a weaker but still significant trend was seen in MERRA2 and VFDv2. While ERA-Interim largely reproduces the trend, WFDEI shows no significant trends. For the considered location, the post-processing of ERA-Interim radiances, based on CRU cloud cover and interannual aerosol loading conducted in the production of WFDEI, has a negative effect on its ability to reproduce the observed trend. The clear-sky type post-processing of ERA-Interim implemented in HySN estimates left the trend close to its original value. The two versions of VFD also differed in their ability to capture the SW↓ trend. This might be due to the maritime location of Bergen and coarser VFDv1 input data for T_min and T_max. A previous study by Parding et al. (2016) showed that circulation type changes could account for a large part of the dimming that was observed in Bergen before around 1980 but a lesser fraction of the subsequent brightening. The fact that ERA-Interim, which does not explicitly account for interannual aerosol changes, picks up the trend while WFDEI, where a correction for interannual aerosol loading has been applied does not, implies that a considerable part of the trend before 2000 must be included in the indirect effects of aerosol changes, which ERA-Interim assimilates, or other factors. On the other hand, MERRA2 accounts for interannual aerosol loading in the time period considered and captures a positive, albeit weaker trend in S_W.

The annual trends in LW↓ during the same period in Bergen were larger in magnitude than those found for SW↓, −8.4 W m⁻² per decade. The observed trend in any calendar month was larger for SW↓, while the LW↓ trend showed more consistency throughout the year. More of the models reproduced the SW↓ trend than the LW↓ trend. Both versions of the VIC type forcing data, VFDv1 and VFDv2, simulated a weak increasing trend in September. Given that the VFD estimates did not produce changes in SW↓ in the same month, the increase is likely due to the clear-sky parametrization and the concurrent simulated increase in September T_d. WFDEI and HySN reproduce the decrease in LW↓ seen in August, while ERA-Interim does not. This points to the fact that changes in near-surface temperature, which are used as a scaling factor and to adjust near-surface humidity in WFDEI and HySN, capture the signal that contributes to the decrease in LW↓. A larger sample of stations measuring incident radiation with a high quality is needed to evaluate how well the models capture trends within the region, particularly given the uncertainty in the observational data.

8 Code and data availability

The HySN data product is available archived in Zenodo (https://doi.org/10.5281/zenodo.1970170, Erlandsen, 2018a). The code used to produce the HySN estimates is written in Python and is available at https://github.com/helene-b-e/HySN.git (last access: 10 June 2019). Further, the particular version of the software code used to produce the HySN estimates validated here is archived in Zenodo (https://doi.org/10.5281/zenodo.1435555, Erlandsen, 2018b). The remaining data sets are available from the various named data providers.

9 Conclusions

Hydrological, ecological, and crop modellers seek landscape-scale data. Norway has a long coastline with mountains, fjords, and small islands. Strong land–sea contrast, high mountains, and a seasonal snow cover that is highly dependent on continentality and altitude results in a fine-scale variability difficult for coarse-scale models to represent. A Python script to downscale and consolidate reanalysis data with high-resolution national gridded temperature data has been developed, which, leaning on previously well-tested empirical relationships, provides estimates of humidity and incident radiation on a fine-scale grid. The downscaled humidity ensures that relative humidity is constrained at 100 %, so that, for example, reasonable evaporation estimates can be sought. The new estimates, HySN, provide humidity estimates with the overall highest quality given for the metrics considered here, also surpassing those based on estimating humidity from temperature alone, such as for the VIC type forcing data. The new estimates outperform the VIC type forcing data and the MERRA2 estimates of incident radiation; however, it is not clear that the new estimates have an added value compared to ERA-Interim and WFDEI. The lack of high-quality historical observations, particularly of incident longwave radiation, hinders a proper evaluation of the data sets.

• Additionally, this study has shown that (ℋ_a) altitude is a significant predictor of humidity, SW↓, and LW↓ in the domain. The coarse-scale estimates of T_d increased on average by 1 ^∘C, SW↓ by 0.7–1.5 W m⁻², and LW↓ increased by as much as 8.6 W m⁻², when adjusted to station altitude.

• Further, the results have shown that a high resolution does not necessarily indicate high-quality estimates. The added value of the high horizontal resolution of the more empirically based estimates does not outweigh the added value of relying on estimates from coarser-resolution numerical weather prediction reanalyses (ℋ_b). Not only is a higher daily temporal correlation (ACC) seen in the estimates based on newer reanalysis data compared to the VIC type forcing data but also a lower mean absolute station bias is seen for several reanalysis-based products (ERA, WFDEI, HySN). VFDv1 and VFDv2 show a 60 % stronger decrease in humidity with distance from the coast than the observations, alluding to the fact that the modified version of the Magnus type formula based on Kimball et al. (1997), implemented in VFD to estimate humidity from daily minimum temperature, is not appropriate for the Norwegian domain. Both VFDv1 and VFDv2 also show a decrease that is several times stronger in solar radiation with latitude than the observations, likely a result of using diurnal temperature range as a proxy for cloud cover, an assumption likely not appropriate in coastal environments and at high latitudes.

To our knowledge reanalysis-based estimates have not been compared with VIC type forcing data for regions within Europe (or Norway specifically). The comparison of model estimates may assist impact modellers that have not yet selected data to use. Some of the findings might help explain persistent errors, for instance found in the timing of snowmelt in a hydrological model. The findings provide emphasis for climate researchers to not only downscale T₂ and precipitation from climate projections, and later use these to estimate humidity and incident radiation, but to utilize the climate model estimates of near-surface humidity and incident radiation. This is already done, for example, in Teklesadik et al. (2017), where the impact of climate change on surface hydrology is examined based on, depending on the hydrological model's structure, bias-corrected climate model output of precipitation, temperature, humidity, and incident radiation from the ISI-MIP project (Hempel et al., 2013). Similarly, we envisage that further work would involve applying the HySN as an input to a hydrological model. Such a model exercise would imply modulating the model's code to accommodate humidity and radiation as input variables. Once the model includes more physically based parameterizations, the sensitivity in simulated runoff to the choice of forcing data can be assessed, including the impact of errors or perturbations in each of the forcing data variables.

The source code for computing HySN has been made available and may easily be configured to use other reanalysis data or other national data sets as input. The compilation of HySN requires merely half a day on a modern desktop computer. Part of the code might also be implemented in a model preprocessor or in the calculation of various indices, so that the variables do not need to be stored for long time spans. Future work entails calculating indices, such as reference evaporation, and updating the input data to ERA5 and a new version of SeNorge once the full historical time series of the two are available. Additionally, sub-daily estimates, the inclusion of terrain features such as slope and aspect, and adding a correction based on the lack of coupling between the land surface and the atmosphere at times when the differences in the local snow cover and snow cover modeled by the reanalysis are large might be promising, as initial results showed that differences in ground snow conditions between the reanalysis and the observations were significant in predicting the difference between ERA-Interim estimates of SW↓ and the observations.