TopoSCALE v . 1 . 0 : downscaling gridded climate data in complex terrain

Simulation of land surface processes is problematic in heterogeneous terrain due to the the high resolution required of model grids to capture strong lateral variability caused by, for example, topography, and the lack of accurate meteorological forcing data at the site or scale it is required. Gridded data products produced by atmospheric models can fill this gap, however, often not at an appropriate spatial resolution to drive land-surface simulations. In this study we describe a method that uses the well-resolved description of the atmospheric column provided by climate models, together with high-resolution digital elevation models (DEMs), to downscale coarse-grid climate variables to a fine-scale subgrid. The main aim of this approach is to provide highresolution driving data for a land-surface model (LSM). The method makes use of an interpolation of pressurelevel data according to topographic height of the subgrid. An elevation and topography correction is used to downscale short-wave radiation. Long-wave radiation is downscaled by deriving a cloud-component of all-sky emissivity at grid level and using downscaled temperature and relative humidity fields to describe variability with elevation. Precipitation is downscaled with a simple non-linear lapse and optionally disaggregated using a climatology approach. We test the method in comparison with unscaled grid-level data and a set of reference methods, against a large evaluation dataset (up to 210 stations per variable) in the Swiss Alps. We demonstrate that the method can be used to derive meteorological inputs in complex terrain, with most significant improvements (with respect to reference methods) seen in variables derived from pressure levels: air temperature, relative humidity, wind speed and incoming long-wave radiation. This method may be of use in improving inputs to numerical simulations in heterogeneous and/or remote terrain, especially when statistical methods are not possible, due to lack of observations (i.e. remote areas or future periods).


Introduction
Simulations of land-surface processes are important for performing assessments of a wide range of earth systems, under current and possible future climates. This task is problematic in complex terrain due to the inter-connected problems of (i) the high resolution required of model grids to capture strong lateral variability caused by, for example, topography, surface or sub-surface processes (e.g. Gubler et al., 2011;Riseborough et al., 2008;Arnold and Rees, 2009); and consequently, (ii) the lack of accurate meteorological forcing data at the site or scale it is required (Thornton et al., 1997;Liston and Elder, 2006). This can be due to the lack of meteorological observations (i.e. spatial coverage, temporal extent/continuity, or variables measured are insufficient for purpose) or lack of representative observations where surface variability is high. Gridded data products produced by atmospheric models can, in part, fill this gap (e.g. Frauenfeld, 2005;Pereira-Cardenal et al., 2011;Akhtar et al., 2009;Vu et al., 2012); however, in many cases not at an appropriate spatial resolution to drive land-surface simulations (i.e. site scale), and therefore require some form of downscaling of variables.
Basic downscaling approaches (also referred to as disaggregation), utilise fixed relationships between the coarsegrid fields of atmospheric models and a subgrid surface, that describes a given variable at higher resolution than the coarse-grid. Such approaches developed from the recognised need to represent subgrid surface heterogeneity in climate models and therefore the requirement to disaggregate the coarse-grid climate forcing to the subgrid land-surface (e.g. Dickinson et al., 1986;Wood, 1992;Koster and Suarez, 1992;Seth et al., 1994). For example, Giorgi et al. (2003) and Dimri (2009) downscale the temperature field according to the subgrid elevation and a fixed lapse rate. In a coupled system it is important that the coarse-grid mean is conserved by such approaches (Giorgi et al., 2003). More complex approaches to downscaling of climate data can be broadly divided into dynamical or statistical methods (Schmidli et al., 2007) which are used to increase the resolution of largescale climate fields (see for full discussion Wilby and Wigley (1997)). Dynamical methods achieve this by using a limited area model at a higher grid resolution e.g. a regional climate model (RCM), to simulate fine-scale processes which are consistent with large-scale climate fields (Giorgi, 2006). While an RCM grid resolution could be increased further, in practice the effective resolution is limited by the complexity of the numerics that must be solved at each model time step, to the order of 10 1 km (Kotlarski and Block, 2005). Statistical methods derive empirical relationships between large-scale predictor fields and local observations (Maraun and Wetterhall, 2010). These methods are computationally efficient but the coverage and effective resolution is often limited by the density of observations, especially in mountainous (i.e. data-poor) areas. Additionally, it is unknown whether empirically derived relationships are valid outside the time window used for calibration. To compliment these approaches there is a growing number of physically inspired, computationally efficient approaches that use physical relationships and high-resolution surface information, that is, digital elevation models (DEMs), to distribute fine-scale forcings (meteorological stations or coarse-grid centre point) over wide areas (e.g. Liston and Elder, 2006;Tarboton and Luce, 1996;Marks et al., 1999;Jarosch et al., 2012). These distributed forcings can then be used for full 2-D, point scale or lumped model simulations.
Another interesting statistical approach presented by Schomburg et al. (2012) utilises empirical relationships between the atmospheric coarse field and high-resolution surface data to disaggregate variables. The main aim of this work was to improve the forcing of soil-vegetation transfer models either in offline simulations or fully coupled model systems. However, the emphasis of this work is on improving the representation of the land surface, via a distributed forcing, within climate models; moreover, the scheme resolution of 400 m would be too coarse to resolve many surface processes in mountain areas.
In complex terrain, topography-based gradients of meteorological variables (i.e. related to elevation, aspect, slope, etc.) can often dominate over horizontal gradients (i.e. latitude/longitude) within a region that is of comparable size to a typical coarse climate cell (e.g. 50-100 km). An example of a method that has successfully encoded this assumption is the PRISM framework (parameter-elevation regressions on independent slopes model) which provides a statistical, topography-based mapping of climate observations (Daly et al., 1994(Daly et al., , 2002. We do not neglect the fact that other forms of heterogeneity may be important in modifying surface fluxes, such as surface cover, but it is important to distinguish between surface heterogeneity that clearly affects atmospheric forcing, e.g. the effect of elevation on air temperature, and those characteristics that become important when performing a coupled land-surface-atmosphere simulation (e.g. the effect of soil properties on near-surface air temperature).
Climate models provide spatially and temporally continuous fields which are physically consistent and therefore are useful tools for forcing regional-scale land-surface studies (Machguth et al., 2009;Kotlarski et al., 2010). In addition, they provide a thorough description of the atmospheric column by providing data fields at many pressure levels between the earth's surface and top of the atmosphere.
The aim of this study is to develop methods that use the well-resolved description of the atmospheric column provided by climate models, together with high-resolution DEMs, to downscale coarse-grid climate variables to a finescale subgrid. The main motivation of this approach is to provide high-resolution driving data for a land-surface model (LSM), understood here as any process-based surface model that simulates the interface of the land surface/subsurfaceatmosphere system. Additionally, the approach should be consistent (methodologically, spatially, temporally) and not reliant on observations. Therefore, the design criteria for this method are (1) it provides high-resolution (< 100 m) downscaling of climate data based primarily on DEM-based information; (2) it is as physically based as possible and therefore has minimal reliance on observations and likely remains valid under future conditions; (3) it employs simple methods which are computationally efficient; (4) it may be used as part of a modelling chain with a lumped representation of the subgrid domain (e.g. Fiddes and Gruber, 2012) for large area applications, as well as 1-D points and 2-D grids. Our approach therefore largely assumes vertical gradients to dominate horizontal gradients within a given model grid box. In this study, we describe this method and its application with ERA-Interim data, a 4-D-VAR reanalysis (3rd generation) which uses the ECMWF climate model, although the method could be equally used with other reanalyses (e.g. NARR, NCEP/NCAR, JRA-55, NASA MERRA, or RCM derived fields). Our methods are then evaluated against a large number of observations over a wide area of complex terrain in the European Alps as well as compared to a set of reference methods. The methods proposed here aim to provide an alternative to statistical methods when observations are not available (remote areas or future periods) and be complimentary to dynamical methods (i.e. they could be used to further downscale RCM output to site scale).

ERA-Interim
ERA-I is a global atmospheric reanalysis produced by the ECMWF. The ERA-Interim data assimilation system contains many improvements both in the forecasting model and analysis methodology relative to ECMWF's previous reanalysis, ERA-40, including the use of 4-dimensional variational analysis, a revised humidity analysis, the use of variational bias correction for satellite data, and other improvements in data handling (Dee et al., 2011). ERA-I provides meteorological data from 1 January 1979 and continues to be extended in near-real time. Gridded products include a large variety of 3-hourly (00:00, 03:00, 06:00, 09:00, 12:00, 15:00, 18:00, and 21:00 UTC) grid-surface fields (GRID) and 6-hourly (00:00, 06:00, 12:00, and 18:00 UTC) upperatmosphere products available on 60 pressure levels (PL) with top of the atmosphere located at 1 mb. ERA-I relies on a 4-D-VAR system which uses observations within the windows of 15:00-03:00 UTC and 03:00-15:00 UTC (in the next day) to initialise forecast simulations starting at 00:00 UTC and 12:00 UTC, respectively. In order to allow sufficient spin-up, the first nine hours of the forecast simulations are not used. All fields used in this study where extracted on the ECMWF reduced Gaussian N128 grid (0.75 • × 0.75 • ). Six PLs are used in this study covering the range of 1000-500 mb (1000,925,850,775,650,500), corresponding to approximately an elevation range of 150-5500 m a.s.l. ERA-I fields used are listed in Table 1.

Pressure levels below model surface
While ERA-I model levels are computed from the model orography surface to top of the atmosphere, pressure-level data is given in the interval 1000 mb (approximately sealevel)-1 mb (top of atmosphere). This means that pressurelevel data exist below the model-grid surface in regions with rugged topography. The extrapolation of fields below the surface uses different methods to those above the model surface. For example, geopotential is extrapolated below the model surface as a function of surface geopotential, surface temperature, temperature at mean sea level, surface pressure and pressure-level value. Temperature is extrapolated below the model surface to a given pressure level as a quadratic function of surface temperature, surface pressure and pressure-level value. Wind and relative humidity are both constant below model surface and equal to the lowest model level values (ECMWF, 2011). Additionally a greater quantity of data is assimilated from above surface observations. Therefore, it can be expected that there is a difference in quality of above/below grid pressure-level data. In addition, measurement locations below grid level are more likely to be in valleys and therefore there will be greater exposure of observations to subgrid processes, not represented by the data. See ECMWF (2011) for full details of extrapolation methods. However, it is difficult to disentangle this effect from the fact that measurements above the model surface are more likely to represent the free atmosphere and, therefore, be more accurately simulated than those strongly effected by topographic effects (e.g. Mesinger et al., 2006;Jarosch et al., 2012) such as inversion layers or topographically modified wind fields. This issue is addressed in the results section.

Evaluation data sets
See Table 1 for an overview (by variable) of the evaluation (OBS) data sets used in this study. The MeteoSwiss automatic meteorological network (ANETZ) covers 40 stations (hourly data) ranging 1132-3580 m a.s.l. and represents both high mountain locations and valleys. The IMIS station network of the Swiss Institute for Snow and Avalanche Research is biased towards high alpine locations (there are few valley stations) but represents topographical heterogeneity in terms of slope and aspect more comprehensively than ANETZ stations. Network elevation range is 1562-3341 m a.s.l. Tenminute measurements from the Alpine Surface Radiation Budget network (ASRB) (Marty et al., 2002) network has 9 stations ranging 370-3580 m a.s.l. All OBS data sources where aggregated to daily mean values to enable comparison with ERA-I fields at a common resolution. An additional analysis on diurnal cycles required aggregation at 3-hourly time steps. See Fig. 1 for the locations all stations used in this study and elevation distribution of stations by variable. See Sect. 6.1 for assimilation issues related to evaluation data sets.

Precipitation climatology
The CRU Alpine precipitation data set (hereafter referred to as CRU) is used as a climatology in the precipitation scheme. It provides monthly precipitation totals, for the period 1800-2003, gridded at 10 arc-min resolution over the Alpine region. The data set is based on 192 long-term homogenized precipitation series from meteorological stations across the study domain and a high-resolution precipitation climatology for the period 1971-1990. Full details are available in Efthymiadis et al. (2006).

Data quality control
OBS values outside acceptable limits were removed automatically by applying physically plausible thresholds to all data sets. Non-changing values beyond prescribed time limits were also screened out (e.g. indicating iced wind propeller). These checks follow the methods of Meek and Hatfield (1994). Discontinuous data sets are valid as the testing strategy follows a point by point comparison between ERA-I Table 1. Pressure-level and surface fields indicate variables which are computed, time steps they are obtained on from ERA-I, sources of validation data (assimilated or non-assimilated), and total stations used in evaluation. Variables in section "other" are used only to support calculations, and therefore columns related to validation are blank. (a) (b) and observations, therefore only data that exists at a given time and location in both ERA-I and observations is carried through to the analysis. This ensures that the maximum possible quantity of valid data was used in the study and makes error-prone gap-filling unnecessary. In aggregation we were careful to ensure that only complete data sets were used in summed values and acceptable levels of datagaps (5 % threshold) were allowed in averaging procedures in the interest of preserving data. No data-gaps were tolerated in summation calculations. Further details are given in the text where relevant.

TopoSCALE
We downscale the variables required to drive an LSM from ERA-I pressure-level (PL) and grid-surface (GRID) fields (Table 1). Input pressure-level fields used are, air temperature (T pl ), relative humidity (Rh pl ), wind components U and V , which are converted to wind speed (Ws pl ). Input gridsurface fields are downwelling global radiation (SW↓ grid ), downwelling long-wave radiation (LW↓ grid ) and precipitation (P grid ). Accumulated values of SW↓ grid and LW↓ grid are converted to time step averages of W m −2 and accumulated P grid is converted to a mean rate of mm h −1 , prior to scaling. The temporal resolution of surface fields is 3 h and PL fields, with native resolution of 6 h, are interpolated to the same 3 h time step (Appendix B). Additionally, the fields T grid , Td grid and SW TOA are used indirectly for radiation computations.
Locations at the coarse-grid level or the fine-scale subgrid are referred to as GRID and SUB, respectively.

Pressure-level fields
Fields derived from pressure levels (T sub , Rh sub , Ws sub ) are computed directly from pressure-level data in two steps: (1) pressure-level elevations (m a.s.l.) are estimated in a standard way by normalising geopotential heights by gravity at sea level (Appendix C1).
(2) Value at SUB elevation is linearly interpolated from data at pressure levels above and below SUB elevation (Fig. 2). Ws sub is derived from U and V wind components after interpolation. Topographically modified wind fields can be additionally computed according to a simple wind sub-model (Liston and Sturm, 1998) which adjusts the speeds and directions according to topographic slope and curvature relationships. To perform the wind modification calculations, the local slope, aspect, and topographic curvature (a measure of relative prominence with respect to surrounding terrain) are required (Appendix C2).

Radiative fluxes
LW↓ sub is computed by deriving a cloud component of allsky emissivity at grid level and using T sub , Rh sub to describe variability with elevation. First, clear-sky emissivity at SUB ( cl sub ) and GRID ( cl grid ) are computed according to Konzelmann et al. (1994): where x1 = 0.43 and x2 = 5.7 (Gubler et al., 2012) and water vapour pressure, pV is a function of Rh (Appendix C3). The all-sky emissivity is computed at GRID using LW↓ grid and the Stefan-Boltzmann equation: as where σ is the Stefan-Boltzmann constant of 5.67 × 10 −8 Js −1 m −2 K −4 . We estimate the cloud-based component of emissivity ( ) at GRID though subtraction of cl grid from as grid in order to apply this correction directly at SUB. Finally, LW ↓ sub can be computed accounting for elevation changes in T and Rh by This approach assumes that cloud emissivity at GRID and SUB elevations are the same, but accounts for reduction of clear-sky emissivity with elevation. This is important as the steepest gradients in LW↓ are often found in clear-sky conditions due to reduction in atmospheric water vapour with elevation. After elevation correction, terrain effects are accounted for by reduction of LW↓ sub by multiplication with the sky-view factor, being the fraction of sky that is visible at SUB (V d ). This assumes that LW↓ is isotropic. SW↓ sub is computed in a three-step process: (1) partitioning of SW↓ grid into direct and diffuse components, (2) elevation adjustment of direct, and (3) topographic correction of both diffuse and direct at point scale. SW↓ grid can be partitioned into direct (SW↓ dir grid ) and and diffuse (SW↓ dif grid ) components according to the hourly regression model of Ruiz-Arias et al. (2010b) which has been developed based on 21 stations in Europe and USA (Appendix ). This method uses the clearness index, which is computed by ratioing SW↓ grid against irradiance at top of atmosphere, SW↓ toa and in doing so estimates a solar transmissivity of the atmospheric column (Fig. 3a). It should be noted that the regression model was developed on hourly data, whereas we apply it to 3 h averages. The vertical gradient of global irradiance between GRID and SUB is mainly determined by the direct component together with difference in the optical path length m, assuming that attenuative properties of the atmosphere are constant between the two elevations. Therefore, we apply an elevation correction to SW↓ dir grid only, largely following methods of Ruiz-Arias et al. (2010b) (Fig. 3b). First, m is computed as where z is difference in elevation and θ z is the solar zenith angle. We can then solve the Beer-Lambert law for direct irradiance to obtain the broadband attenuation coefficient (k): where m = 1/ cos θ z (except for large values of θ z ). The difference in SW↓ dir due to elevation difference between GRID and SUB can then be found as Equation (6) shows direct irradiance should increase exponentially with elevation given constant k and elevation-based change in irradiance is maximum when sun is at zenith and zero when sun is at horizon. As the correction is only applicable to clear-sky conditions, it is applied when the air-masscorrected clearness index k t is greater than 0.65 (Perez et al., 1990).  (3) components through a clearness index by ratioing against extraterrestrial radiation (1). (B) An elevation correction is applied to SW↓ dir grid (1) to obtain SW↓ dir sub (2) based on z , M and θ z . (C) Topographic correction is applied accounting for illumination angle, θ i , horizon angle, θ h and sky view factor V d .
SW↓ dir sub is computed by correcting for terrain effects of slope, aspect and horizon at SUB according to Dozier and Frew (1990) and Dubayah and Rich (1995) (Fig. 3c). This is achieved by computing the illumination angle i, which is the angle the incident direct beam makes with the slope normal and varies with solar zenith θ 0 and azimuth angles φ 0 , and local slope angle S and aspect A: By ignoring variation in latitude and longitude within a given grid box θ z and φ 0 can remain constant, which over short-length scales is a reasonable simplification (e.g. Dozier and Frew, 1990). As slope = 0 at GRID, cosi grid is simply given by Additionally, cast-shadows and self-shadowing effects are often important in complex terrain and are accounted for through local horizon elevations. Wherever cosi sub is negative, the point is self-shadowed, that is, the sun is below the horizon formed by the local slope and SW↓ dir sub is set to 0. Cast shadows are found by horizon elevations and given as δ, a binary shadow mask. Topographically corrected SW↓ dir sub is then given by first removing the GRID cosine correction and then multiplying by SUB cosine correction: where horizon elevations are either explicitly given in n directions for 1-D/2-D simulations or parameterised for use with a lumped scheme (e.g. Fiddes and Gruber, 2012) as a function of local slope and V d in order to detect shading. Computation of SW ↓ dif sub assuming isotropy requires only V d :

Precipitation
Precipitation patterns in mountain regions are driven by a range of complex mechanisms which depend upon, for example, season, geographical climate (maritime, continental) and structure of orography (Leung and Ghan, 1995;Lundquist et al., 2010;Smith and Barstad, 2004). Precipitation is therefore strongly variable in both time and space. Figure 4 gives an overview of the combined climatology and lapse rate approach used to address this problem (e.g. Früh et al., 2006). It should be noted that this routine can also be implemented using only the lapse rate, in order to remove any dependency on climatology data. We acknowledge that the quality of the data set is heavily dependent on density of observations and therefore, better than average results should be expected in our study domain. First, the elevation signal of the climatology data is removed by dividing by the non-linear lapse rate (λp) of Liston and Elder (2006) (Appendix C5), and then normalising to the GRID reference elevation. P grid is then disaggregated according to the subgrid variability as described by the normalised climatology (now elevation independent). Each of the ith (in this study 1-25) climatology grid-cells P clim are normalised by the sum of P clim contained in each ERA course grid box, resulting in a subgrid disaggregation factor W sub , In this study we use the CRU climatology but other sources of subgrid observations could be used. Globally, or near globally available data sets (gauge and satellite based) include the Tropical Rainfall Measuring Mission (Huffman et al., 2007), Global Precipitation Climatology Centre (Beck et al., 2005) and Climatic Research Unit (New et al., 2002) or regionally available such as PRISM (Daly et al., 1994) or even direct observations if available. The product of P grid and W sub generates a subgrid distribution of precipitation at the climatology resolution, that is conservative of the coarsegrid forcing P grid . Finally, λp is applied to capture fine-scale precipitation-elevation gradients in order to obtain P sub , It should be noted that this method simply applies a scaling to the precipitation field and does not estimate the precipitation phase (snow-rain partition), this is done by the model which the precipitation field is used to drive.

Reference methods
The following is a description of the reference parameterisations with which we compare the current scheme. We do not intend this comparison to be exhaustive, but to merely serve as a reference point, based on common parameterisations.
T grid is simply extrapolated according to a fixed lapse of 6.5 • C km −1 (e.g. Blandford et al., 2008). Rh is not a linear function of elevation and so the relatively linear dew point temperature (Td) is often used (Liston and Elder, 2006). Rh can be converted to Td as a function of T and pV. In this study Td grid is available, so this step is unnecessary.
Td sub can be computed using a variable lapse rate (Kunkel, 1989): where λ is a vapour pressure coefficient that varies during each month of the year (Kunkel, 1989) and constants b = 22.452, c = 272.55 • C are given by Buck (1981). Finally, Td sub is then converted back to Rh sub as a function of T sub . LW↓ is parameterised as a function of T and pV and cloud cover according to the clear-sky formula of Konzelmann et al. (1994) (Eq. 1) and the all-sky formula of Pirazzini et al. (2000) while accounting for V d , where cl is given by Eq. (1), N is given by the ERA-I total cloud product (0-1), p1 = 6, p2 = 4 and as = 0.979. Wind is adjusted by 40 % per km (i.e increased above grid level reduced below grid level) (Plüss, 1997). SW↓ is not compared to a reference method as the methods (i.e. partitioning, elevation and topographic correction) used in this study are commonly used elsewhere (i.e. Oliphant, 2003;Schroeder et al., 2009). Precipitation is scaled with a non-linear lapse rate (Liston and Elder, 2006).

Location
The study region contains the entire Swiss Alps and uses 19 ERA-I grid boxes to cover all observations that are used to evaluate the methods. Switzerland contains one of the most densely observed mountain regions in the world and, therefore, is a suitable region within which to evaluate this method and specifically its performance with respect to vertical information, as it covers a large elevation gradient of 195-4634 m a.s.l. (Fig. 1).

Set-up
The experimental strategy is as follows: results of the current methods (TopoSCALE) are compared to (1) unscaled grid level ERA-I fields (GRID) and (2) 2.3). Statistical evaluation is primarily performed using the correlation coefficient (R), the root mean squared error (RMSE) and bias (BIAS) in order to test for systematic errors, expressed simply as BIAS = (sim − obs).

Results
In this section results are presented as follows: (1) pressurelevel-based results (T , Rh, Ws, LW↓); (2) surface-based results (SW↓ and precipitation); (3) seasonal error signatures; (4) diurnal cycles; and (5) elevation effects, both absolute and relative to grid level. For Rh (210 stations), TopoSCALE gives a significant improvement in the correlation and modest improvement in RMSE (due to the poorly performing cluster at high humidity). There appears to be a high degree of uncertainty in saturated conditions (i.e. measurements at or close to 100 %) in all cases. TopoSCALE shows significant improvement over GRID and REF (approximately the same performance) particular in the interval 0-60 % which is significant for processes such as sublimation which occur in dry atmospheric conditions. Both GRID and REF seem unable to represent humidities less than 30 %. This could possibly be that dry (often well below 30% humidity) Foehn winds are represented in the pressure-level data but absent from surface data.

Pressure-level-based results
Several discontinuities were observed in the wind time series that are possibly related to changing patterns of data assimilation (Fig. 7a). At least one of these artefacts fits to a major data introduction (European wind profilers in 2002). Therefore, the wind analysis was restricted to a three year period (1996)(1997)(1998) that was stable. Comparison of distributions of wind speed show a large improvement of TopoSCALE over GRID especially in the R value and BIAS. There is still a large degree of scatter especially at high wind speeds. The most significant result is an improvement in BIAS with a shift of the distribution to the 1 : 1 line, while the GRID data is not able to represent values greater than 5 m s −1 . Additionally, Fig. 6 gives a comparison of TopoSCALE and TopoSCALE + wind sub-model for the ANETZ station at Natschen above Andermatt. The 2-pronged error signature (i.e. both topographic wind speed reduction and enhancement) is corrected by the sub-model. The wind sub-model is difficult to test widely (all other wind results do not include the sub-model) as topographic location data is often not precise enough for point validation, especially where locations are peaks or ridges (i.e. flat ridge can be extracted as a 45 • north face even on a 25 m DEM, e.g. Fisher et al., 2004). This station was chosen for its position on a large slope well represented at the DEM resolution.
For LW↓ (9 stations Figure 5 gives the density scatter plot for SW↓ grid OBS (9 stations) against GRID and TopoSCALE results (MOD) for both all-sky conditions and clear-sky conditions (defined as k t > 0.6). Best RMSE performance is seen in clear-sky conditions due to removal of a large proportion of cloud-based uncertainty. However, BIAS is higher in clear-sky conditions due to residual elevation effects affecting the larger direct beam component. TopoSCALE reduces this BIAS by around 3.5 W m −2 as well as improving the RMSE score by  1999-2001 and 2005-2009. Therefore, the wind analysis was restricted to a three year period (1996)(1997)(1998)) that was stable. (b) A common problem with climate models is a low number of dry days which is compensated for by too high frequency and low intensity precipitation. The percentage dry days of OBS is much higher and distribution not even overlapping that of GRID data (IMIS-Data © 2013, SLF). 5 W m −2 . BIAS is also reduced under all-sky conditions by TopoSCALE but more modestly and RMSE score is roughly equal. This indicates that TopoSCALE improves the direct beam component most as corrections focus on this part of the radiation budget. In a separate analysis the Erbs partitioning scheme for direct and diffuse SW↓ was tested. As expected the partitioning scheme adds some uncertainty (i.e. R reduced from 0.88 to 0.81/0.75under all-sky conditions). Results for direct/diffuse are negatively/positively biased (−20.4/17.2 %). Full details are not given here as this topic is well covered in the literature (e.g. Erbs et al., 1982;Ruiz-Arias et al., 2010a). Despite high uncertainty introduced in decomposition models, the reaggregation of solar components after elevation/terrain correction minimizes the potential effects in the final terrain corrected estimates (Ruiz-Arias et al., 2010a). Two quantities are important in modelling precipitation, quantity and frequency, the former representing total inputs and the second controlling distribution of those inputs over a given period of time. Figure 7 shows a common problem with climate model precipitation fields -that is "constant drizzle" (i.e. a low number of dry days which is compensated for by high frequency and low intensity precipitation) (Piani et al., 2009;Manders et al., 2012). The percentage dry days of OBS is much higher (and not even overlapping) that of ERA-I GRID data. This cannot be changed by the current scheme as precipitation can not be created (to be conservative), but only distributed to the subgrid according to the scheme. A conservative approach would require a temporal redistribution of precipitation as opposed to the spatial corrections we apply in this study.  Figure 9 also highlights the improved simulation of both extremes.

Diurnal cycles
The diurnal cycle in surface and boundary-layer variables is important for the global climate system (Dai and Trenberth, 2004), and particularly in simulating daily variation in the surface energy balance. Figure 10 shows the diurnal cycle of SW↓ and T , two fields characterised by distinctive diurnal cycles, in order to investigate the performance of the scheme at sub-daily timescales. Additionally these fields represent surface (SW↓) and pressure-level (T ) fields. We calculated the average of all 03:00-00:00 UTC 3 h time steps over the entire study period for months of December and June. A subset of OBS stations is presented, representing an elevation range of 370-3580 m a.s.l. In general, the diurnal cycle of SW↓ appears to be well reproduced by ERA-I. However, seasonal differences are apparent with more accurate simulation in June than December as indicated by the full range of values at 12:00 UTC being more comprehensively represented. Lower amplitudes of diurnal cycles in T make the analysis less clear. In December, diurnal cycles are quite strongly smoothed at low elevation whereas cycles are virtually non-existent in the OBS data at high elevation. In June there is a degree of smoothing at low elevations but cycles are generally reproduced. However, at high elevation there is a very strong smoothing. Where TopoSCALE performs less well for T (i.e. winter and high elevation) is likely related to poor representation of surface boundary layer in ERA-I data (cf. Sect. 6.3). Figure 11 gives boxplots of deviation of daily values of MOD from OBS, as defined above and grouped according to month of the year in order to investigate seasonal patterns in the error signature.  Figure 12 shows the daily mean error of GRID, REF and TopoSCALE results (MOD) with respect to OBS as a function of relative elevation of station (e.g. station elevationgrid elevation). Each box may contain multiple stations as long as they share the same elevation difference from their respective ERA-I grid cell. The plot is binned into 400 m intervals (300 m for radiation). This analysis was performed to investigate any elevation dependency of the error signal as well as to look at the effect of the different methods implemented in the ERA-I model to compute variables on pressure levels above and below the model surface (cf. Sect. 2.2). The red box represents surface data (grid-level ±200 m) in order to investigate the relative performance gain/loss close to the grid-surface. This may point to suitability of TopoSCALE outside of mountain areas.

Elevation effects
The results for T show larger error for stations below grid (RMSE = 2.46) than above (RMSE = 1.60). This result is also slightly negatively biased. This shows that the extrapolation of data below grid level produces a poorer result and only slightly better than REF. The fact that observations tend to be colder indicates that non-represented subgrid effects, such as inversions, could be significant in driving this bias. Above grid level there are large improvements over REF and GRID. REF shows the expected result that error related to lapse-rate-based approaches increases with the distance over which they are applied. The Rh plot shows that TopoSCALE is increasing positively/negatively biased with distance above/below grid level as compared to REF.
However the absolute magnitude of error is much lower. Ws bias in GRID error signature is corrected by TopoSCALE (albeit slightly overcompensated above grid). TopoSCALE performs better below grid level, possibly due to lower absolute Ws magnitudes leading to lower error values.
Looking only at surface data (red-box), pressure-levelbased results (i.e. TopoSCALE) outperform surface databased results (i.e. GRID and REF) in all cases but most significantly in T and Ws. This is quite surprising as it would seem likely that the surface data should contain more of the boundary layer effect (cf. Sect. 6.3). However, this indicates that TopoSCALE could also be usefully applied in locations close to grid level without reduced quality over surface-based data.

Discussion
Reanalyses are complex products in that they combine a climate model with observations. This section provides a discussion of key issues relevant to the use of reanalysis and other climate data sets, in order to place this study and results in context, as well as to highlight some important limitations to this approach.

Assimilation issues
Reanalyses assimilate a large number of observations in spatially and temporally varying quantities and densities. It is therefore important to know which observational data sets are assimilated into the ERA-I product, as this not only affects how independent observations are in terms of validation, but will also suggest how the performance of ERA-I (and therefore TopoSCALE) can vary with observation density. Assimilated data that is also used for evaluation in this study originates from the SYNOP registered MeteoSchweiz stations (a subset of the ANETZ network), and only affects observations of air temperature and relative humidity (Table 1). Furthermore, for screen-level analysis (2 m temperature, 2 m relative humidity) surface observations that differ by more that 300 m from the model orography are rejected in the ERA-I assimilation.

Bias and spatial-temporal variability of errors
LSM results are sensitive to bias in climate data (e.g. Berg, 2003), which may result in unrealistic estimates of mass,  Fig. 9. Performance of TopoSCALE precipitation at monthly and annual scales compared to GRID and REF. The description of the spatial distribution of precipitation included in TopoSCALE gives improvements over a purely lapse-rate-based approach (REF). energy, momentum exchanges between the atmosphere and surface (Maurer et al., 2001a, b). Therefore, bias correction is usually regarded as a crucial step in providing accurate driving fields to a land-surface or impact model (Hagemann et al., 2011). In this study we have chosen to conceptually separate "bias" and "scale" (often combined in downscaling routines based on station data) in order to focus on the problem of topography-based scaling, as this is not reliant on observational data sets. Therefore, the treatment of bias can be performed in a second step with a reduced influence of scale differences. However, we acknowledge that bias correction is often necessary to provide accurate fields to surface models. Reanalysis can be seen as an imperfect model combined with incomplete data and output should not be equated with "observations" or "reality". The changing mix of observations, and biases in observations and models, can introduce spurious variability and trends into reanalysis output. Observational constraints, and therefore reanalysis reliability, can vary considerably depending on the location, time-period, and variable considered. Another problem is that mixing observations with models tends to violate conservation laws. Most significant to this study is the fact that reanalysis will likely be closer to reality at locations with higher observation densities (i.e Europe and specifically the European Alps, in contrast to other high-mountain regions).

Subgrid issues and boundary-layer effects
Due to the coarse resolution of current reanalysis data sets (typically 0.75 • -1.5 • ), various processes are unresolved by the model. An important example is temperature inversion in mountain valleys, which will not be captured by the data. Another important process is local-scale rain shadows caused by unresolved topographic barriers and shallow convection which is parameterised by a bulk mass flux scheme in the ECWMF model (Tiedtke, 1989), and therefore cannot resolve the level of spatial differentiation that is present in surface measurements. The surface boundary layer (as opposed atmospheric boundary layer) will have a residual effect upon surface measurements, which will not necessarily be present in pressure levels representing the free atmosphere. For example, turbulent exchanges of sensible heat fluxes can be a significant contributor to energy exchange between surface and atmosphere (Cline, 1997;Helgason and Pomeroy, 2012). These effects will also likely affect diurnal cycles of observations. However, the magnitude of these effects is not quantified by this study.

Conclusions
This study has proposed a method that can efficiently provide meteorological variables to an LSM operating at high resolution in complex terrain. In addition, it provides a means to generate driving data in remote areas due to non-reliance on measurements. The schemes focus is on variables that can be derived from pressure-level data, however, surface fields are also computed in order to provide a consistent set of driving meteorology required by an LSM. Important limitations of the approach are described in Sect. 6 but can be summarised as related to (1) assimilation issues (i.e. possible reduction in performance in data-poor areas), (2) reduced performance below grid-level (although this is not universal in gridded data sets), (3) bias in gridded climate data, and (4) subgrid phenomena that are not resolved by the input data (such as temperature inversions). Specifically, in terms of variables computed in this study, strong improvements in the subgrid radiation scheme would likely result from the availability of radiative fluxes (LW↓, SW↓ dir and SW↓ dif ) on pressure levels. The benefit of such an improved description of vertical profiles and diffuse/direct partitioning would be of great relevance as large areas globally are subject to rugged topography (cf. Gruber, 2012;Körner et al., 2011;Meybeck et al., 2001). As an outlook, partitioned SW↓ components are now archived in the current ECWMF operational model (ECWMF, personal communication, 2012), so will likely be available in future generations of reanalysis, but possibly only at the surface. Precipitation could be improved by a more rigorous sub-model such as that proposed by, for example, Smith and Barstad (2004). We attempt to account for subgrid orographic effects such as rain-shadows by implementing a description of variability through a climatology, but this is dependent on the quality of the climatology data set. However, the next generation of reanalysis (e.g. ERA-20C) are likely to deliver large improvements in this respect due to higher model resolutions improving the representation of orographic precipitation. Additionally, when used with RCM (e.g. CORDEX, the first globally integrated RCM project) or weather model it is expected that precipitation-based error would be reduced significantly. In sum, the core strengths of the scheme we have described in this study are the following: 1. Generally demonstrates improved skill in downscaling climate data in complex terrain, as compared to reference methods and good performance when evaluated against measurements. This result is most pronounced in pressure-level variables.
2. Provides a means to generate downscaled data when statistical methods are not possible i.e. in remote, datapoor areas or future time-periods.
3. Provides spatially and temporally continuous meteorological fields at point-scale which are physically consistent.
4. Is efficient and therefore can be used to derive long time series or data over large areas.

C5 Precipitation lapse rate
The non-linear lapse rate (λ p ) used to calculate precipitation in REF and TopoSCALE after Liston and Elder (2006): where precipitation factor, pf = 0.27 (mean of monthly values given by Liston and Elder (2006)) and elevation difference between GRID and SUB is given by eD. Elevation adjusted P sub is then computed from (disaggregated) P grid :