Modelling forest snow processes with a new version of WaSiM

ABSTRACT We present a new model extension for the Water balance Simulation Model, WaSiM, which features (i) snow interception and (ii) modified meteorological conditions under coniferous forest canopies, complementing recently developed model extensions for particular mountain hydrological processes. Two study areas in Austria and Germany are considered in this study. To supplement and constrain the modelling experiments with on-site observations, a network of terrestrial time-lapse cameras was set up in one of these catchments. The spatiotemporal patterns of snow depth inside the forest and at the adjacent open field sites were recorded along with snow interception dynamics. Comparison of observed and modelled snow cover and canopy interception indicates that the new version of WaSiM reliably reconstructs the variability of snow accumulation for both the forest and the open field. The Nash-Sutcliffe efficiency computed for selected runoff events in spring increases from −0.68 to 0.71 and 0.21 to 0.87, respectively.


Introduction
Forest cover is a major cause of the spatial variability of snow accumulation and snowmelt, especially in subarctic environments and subalpine mountain regions with a high proportion of forest areas. Forests can thereby store considerable amounts of snow. Pomeroy et al. (1998) reported that boreal forest canopies can intercept up to 60% of cumulative seasonal snowfall. Published values for the maximum interception storage capacity of a conifer forest range from a few millimetres (Hedstrom and Pomeroy 1998) to 40 mm of snow water equivalent (SWE) (Storck et al. 2002), depending on snow and climate conditions. Generally, trees have a higher interception capacity for snow than for rain (Lundberg and Halldin 2001). Intercepted snow may be subject to high sublimation rates depending on atmospheric conditions. The amount of sublimation loss depends on the exposure time of snow in the canopy and can rise up to 30% of cumulated snowfall. The remaining 70% of snow in the canopy can be released from the branches as meltwater drip or mass unload , leading to highly non-uniform snow depth patterns underneath the canopy (Buttle and McDonnell 1987, Hedstrom and Pomeroy 1998, Gelfan et al. 2004, Floyd and Weiler 2008, with generally reduced amounts of snow in the forest (e.g. Lundberg et al. 1998, Pomeroy et al. 2002, Winkler et al. 2005, Strasser et al. 2011.
Numerous studies have shown the influence of a forest canopy on the snowmelt processes on the ground (e.g. Faria et al. 2000, Pomeroy et al. 2003, Gelfan et al. 2004, Niu and Yang 2004, Liston and Elder 2006, Musselmann et al. 2008, Ellis et al. 2010, Garvelmann et al. 2014, Gouttevin et al. 2015. Both delaying and accelerating effects of snowmelt have been described (Link and Marks 1999a, 1999b, Tribbeck et al. 2004, Strasser et al. 2011. Decreased global radiation due to shading by trees and reduced wind speed inside the forest may delay snowmelt, whereas increased longwave emission by the trees accelerates snowmelt and becomes an important component of the surface energy balance (Link et al. 2004, Sicart et al. 2004, Essery et al. 2008, Pomeroy et al. 2009). Which opposing effect dominates depends on several factors, such as canopy density, gap size and distribution, geographical location, aspect and meteorological conditions (Pomeroy et al. 2002, Jost et al. 2007, Strasser et al. 2011. For the meteorological variables air temperature, relative humidity, wind speed and incoming shortwave radiation, a considerable heterogeneity can be observed beneath the trees (Withaker and Sugiyama 2005). The simulation of melt rates and the timing of meltwater release from forested areas can hence be significantly improved by including local meteorological observations (Marsh et al. 2003, Strasser andEtchevers 2005). Since these are mostly not available, transfer algorithms are widely used to calculate sub-canopy climatic conditions from open-field site measurements (Strasser et al. 2011). A precise representation of sub-canopy climatic conditions is particularly important for the simulation of the snow surface energy balance at the forest ground.
The effects of the forest canopy on snow accumulation and snowmelt processes need to be considered in simulations of the hydrological response of catchments with a significant fraction of forested area. This holds true especially for the subalpine elevation zone, where the winter season is typically characterized by multiple snow accumulation and snowmelt periods, as is the case in our test sites. These conditions, i.e. a frequent change of the dominant processes, are more challenging to simulate than snow processes in colder regions, such as arctic or high alpine environments, where usually only one long accumulation period and one single melt period occur during winter (Rutter et al. 2009). During the experimental phase of the SnowMIP2 study (Rutter et al. 2009), the simulation of surface energy exchange and snow dynamics in forested environments in particular was investigated. The results showed the need for improved models for forested environments to simulate separately the snow cover accumulation and ablation inside and beneath the canopy.
For the study presented here, we applied the physically-based hydrological model WaSiM (Schulla 1997(Schulla , 2015 to simulate the water fluxes from the atmosphere to the vegetation, further to the soil and finally along the river to the gauge at the outlet of the catchment (www.wasim.ch). WaSiM has been used successfully in a wide range of applications (e.g. Schulla 1997, Niehoff et al. 2002, Kunstmann and Stadler 2005, Cullmann et al. 2006, Wriedt and Rode 2006, Bormann et al. 2007, Kraller et al. 2012, Natkhin et al. 2012, Warscher et al. 2013. The hydrological processes are computed on a regular grid and in hourly or daily time steps. Most recently, WaSiM has been extended for processes typical in specific snow processes, e.g. lateral snow transport (Warscher et al. 2013), and a layered snow model and heat transport in the soil and the snowpack for representing permafrost (Daanen and Nieber 2009). The processes at the snow surface can be computed by means of the energy balance approach. However, snow-canopy interactions, including snow interception processes and the modification of the meteorological conditions inside the forest, are still missing in the currently available version of the model (9.09.01). In catchment-scale hydrological modelling, this type of snow-canopy interaction simulation has rarely been applied so far (one remarkable example is presented in Gouttevin et al. 2015).
In the present study, we describe the implementation of a snow-canopy interaction model to simulate the climatic conditions inside the forest and the processes of interception, sublimation and melt unload of snow on the trees. The used snow-canopy interaction model is based on the work of Liston and Elder (2006) and has already been integrated in a few snow and hydrological models, differing in some of the parameterizations used and level of detail: Strasser et al. (2011) integrated the snow-canopy interaction model into the hydroclimatological model AMUNDSEN to simulate the effect of different types of forests on the dynamics of the mountain snow cover; Förster et al. (2014) transferred the snow-canopy process descriptions into the semi-distributed hydrological model Panta Rhei, and applied it for a forested mountain catchment in the Harz Mountains (Germany); and Marke et al. (2016) made it available as a spreadsheet-based snow model for the (point) location of a meteorological station. Now, we have added the snow-canopy interaction model as a new module to the fully distributed modelling system WaSiM whichother than the mentioned modelscombines both a distributed representation of hydrological processes and detailed descriptions of the processes in the unsaturated and the saturated zones for runoff generation and channel streamflow. The simulations are supported with observations collected by means of a terrestrial time-lapse camera network. Since WaSiM is an open source model, the new model extension will be made freely available.
The main objectives of our study are: • To present an extended version of the water balance model WaSiM including a detailed description of the new model algorithms for snow interception and climatic conditions inside the canopy. • To evaluate the performance of the snow cover simulation in a forested subalpine catchment by comparing modelled state variables with insideforest observations recorded using a terrestrial time-lapse camera network. • To validate the model performance in simulating basin streamflow in two catchments, and to quantify the effect of the snow-canopy extension on the modelled discharge.

Study sites
The study was carried out in two study areas: the Brixenbach valley, a small subalpine catchment situated in the Kitzbühel Alps in Northern Tyrol, Austria (Fig. 1); and the Sieber catchment in the Harz Mountains, Lower Saxony, Germany (Fig. 2). The size of the Brixenbach catchment is 9.3 km 2 , with a mean elevation of 1370 m a.s.l. The highest point (Gampenkogel) has an elevation of 1956 m a.s.l. and the discharge gauge of the Hydrographic Service of Tyrol (installed in 2004) at the catchment outlet is at 818 m a. s.l. The mean annual precipitation sum at the precipitation gauge at Nachtsöllberg (990 m a.s.l.) close to the catchment outlet is about 1400 mm, and the mean duration of snow cover amounts to 132 days (1990-2010. The bedrock belongs to the Paleozoic Greywacke zone and is thus dominated by porphyroids and shales (slightly metamorphic sand-, siltand claystones), partly overlain by Mesozoic dolomites. Mostly shallow cambisols, podsols, partly gleysols andin dolomite areasrendzinas have developed on the Quarternary sediment coverage (moraines, talus deposits, colluvium). The catchment area is mainly covered by oligotrophic cattle pastures (44%) and forests (35%). Rock faces and talus slopes cover 14% of the catchment, and only small areas are used as hay meadows for settlements, ski-slopes and forest roads (Meißl et al. 2017). The forests are dominated by conifers, with spruce (Picea abies) being the predominant tree species. Larch trees (Larix), firs (Abies), mountain pines (Pinus mugo), Swiss stone pines (P. cembra), grey and green alders (Alnus incana, A. viridis) occur in smaller proportions.
The Sieber catchment is located in the Harz Mountains, a low mountain range in the northern part of Germany. The catchment upstream of the Pionierbrücke gauging station (338 m a.s.l.) covers an area of 44.4 km 2 . The Bruchberg in the northern part of the catchment is the highest mountain (927 m a.s.l.). Above the quartzite, greywacke and granite bedrock, sandy loam and loamy sand soils prevail. Besides deciduous forest, meadows and upland peat bogs, the major part of the catchment is covered by coniferous forests, in which Norway spruce (Picea abies) is the predominant wood species (Förster et al. 2014). A complex system of channelsthe Upper Harz Water Management Systemwas constructed in medieval times in the course of mining activities. This is relevant from a hydrological perspective, since water is redirected via the channels across the watershed boundary. Long-term recordings at the Sieber rainfall station (see Table 1) amount, on average, to 1270 mm year −1 . Observations obtained from non-recording raingauges (totalisators) at higher elevations in the Sieber catchment even suggest higher annual totals of more than 1500 mm (Förster 2013). Even though the elevation of the Harz Mountains is low compared to the Alps, the accumulation of a seasonal snow cover is typical in the higher-elevation bands (Förster 2013).

Input and validation data
As meteorological input data for the WaSiM simulations in the Brixenbach catchment (precipitation, air temperature, relative humidity, global radiation and wind speed) the INCA dataset of the Austrian Meteorological Service (Zentralanstalt of Meteorologie and Geodynamik, ZAMG) was used for the period 2009-2016 (Haiden et al. 2011). These data have a temporal and spatial resolution of 1 hour and 1 km 2 , respectively. The first modelling results revealed that convective precipitation events are not always appropriately represented in these data; hence, for precipitation, we replaced the INCA fields with the precipitation recordings of ZAMG and HD Tirol, who operate quite a dense network in the vicinity of the Brixenbach valley (see Table 1 and Fig. 1, small map on the left). Interpolation was done by the inverse distance weighting method with elevation dependence (Schulla 1997). In particular, the main improvement in model performance was achieved by considering the raingauge at Talkaser Alm, situated inside the catchment border (see Table 1).
Similar gridded meteorological data were not available in the Sieber catchment. Instead, hourly observations of air temperature, relative humidity, wind speed and global radiation from the Braunlage automatic weather station (607 m a.s.l.) were used as forcing data. Precipitation was interpolated in a similar way using station data (Table 1 and Fig. 2).
For the Brixenbach catchment, spatial data used for the WaSiM parameterization were resampled to a   (Copernicus 2017), with slight adaptations (e.g. natural grasslands were reclassified as pastures). The digital terrain model ALS-DTM provided by the Federal Province of Tyrol was produced on the basis of Airborne Laser Scanning data from 2008/09. Similar spatial data were also available for the Sieber catchment (Müller et al. 2017). In contrast to the Brixenbach catchment, a spatial resolution of 150 m × 150 m was used due to the larger catchment size. For model calibration and validation, we used the recordings of the discharge gauges at Brixen im Thale (gauge no. 202663, 819 m a.s.l., operated by HD Tirol 2004-2016) and Pionierbrücke (gauge no. 4882161, 338 m a.s.l., operated by Harzwasserwerke). In order to validate the new snow-canopy parameterizations, a network of terrestrial time-lapse cameras along with snow poles was deployed in the Brixenbach catchment to continuously monitor snow depth and interception dynamics (Garvelmann et al. 2013). In the winter 2015/ 16, six cameras were installed in the study catchment at different elevations (Table 2), with pairs of cameras close to each other, one located underneath the forest canopy and the other in the adjacent open field. The camera type used was a standard waterproof outdoor camera (Dörr SnapShot Mobil 5.1) that was originally developed for wildlife observations. The cameras were operated independently of any additional power source in all locations. At each camera site, a wooden snow pole was installed in the camera field of view, with alternating, 10-cm black and red bars as snow depth indicators (as shown in the "Results and discussion" section, Fig. 5). The cameras were set to take one digital image every hour. For the analysis, one picture per day (around noon) was selected to estimate daily snow depth using a semi-automatic procedure described in Garvelmann et al. (2013). The snow interception evolution in the forest canopy was determined from the digital images by qualitative analysis, applying a semi-automatic routine with a specifically developed image processing software that allows one to decide whether intercepted snow is present or not. Furthermore, the pictures from the time-lapse cameras provide useful information about the weather conditions in the area, the snow-covered area in the camera field of view, as well as the state of precipitation.
Weekly snow density recordings were available from the Hydrographic Service of Tyrol for Kelchsau, Jochberg and Kössen, approximately 5-30 km from the study catchment. The snow density values were linearly interpolated to daily values and used to calculate the snow water equivalent (SWE) from the continuous snow depth observations at the camera locations. We assumed this approach to be valid for our application, since snow density can typically be associated with a much lower degree of spatial variability compared to snow depth.

Model description WaSiM
The WaSiM model is a deterministic, physically-based hydrological model that has been designed to compute hydrological processes on a regular grid. The first applications of the model focused on studying landuse and climate change effects for mesoscale catchments (e.g. Schulla 1997, Niehoff et al. 2002. While early releases of the model were capable of simulating natural catchments, recent progress in its development have resulted in a plethora of new features regarding effects subject to water management strategies, e.g. irrigation or reservoir operation. The simulation of processes such as silting up, which considers coating and crusting of the upper soil or glacier evolution, has been refined. The concept of the model allows the selection of different representations for the same process and thus enables the testing and comparison of multiple variants of the model for one application. Different sub-models can be selected by the user for several processes, depending on the availability of data and requirements of the specific application (Schulla 2015). The model version used in this study also includes the recent model extensions relevant for applications in snow hydrology, i.e. computation of snowmelt by means of the energy balance (Warscher et al. 2013), and permafrost (Daanen and Nieber 2009).
A full list of available process descriptions can be found in Schulla (2015). Here, only the most relevant process descriptions used in this study are further explained. A threshold of 0°C with a transition air temperature of ±0.5°C for mixed precipitation is used in WaSiM for the phase change between rain and snow. This was applied in the model set-up for both study catchments. Rain interception is computed using a bucket approach, capable of representing layered vegetation. This approach allows interception in the canopy, and subsequent interception of water in the understory below the trees to be computed. Intercepted water is subject to evaporation. The potential evapotranspiration is computed using the Penman-Monteith approach (Monteith 1965). Infiltration is computed according to Green and Ampt (1911) and is coupled with a complex numerical scheme used for computing the soil water balance based on the Richards (1931) equation. The soil is represented by a one-dimensional discretization of soil layers. The runoff components are computed by balancing the soil moisture water content with infiltration, deep percolation, actual evapotranspiration and saturation. Surface runoff is the water not infiltrating the soil, and interflow is a fast runoff component originating from the unsaturated zone. Even though a two-dimensional groundwater model can be coupled to the numerical soil water scheme, a linear reservoir approach for computing baseflow can be used alternatively (as in this study), representing the outflow from saturated layers. Finally, runoff routing is based on the Manning-Strickler equation (Strickler 1923).

Beneath canopy climate
The new model extension also includes a modification of the sub-canopy meteorological conditions since they differ significantly from those in the open. Beneath the trees, the shortwave radiation, precipitation and wind speed are reduced, longwave radiation and humidity are increased, and the diurnal course of air temperature is attenuated (Link and Marks 1999a, 1999b, Tribbeck et al. 2004, Strasser and Etchevers 2005. As measurements of the meteorological variables are usually taken at observation sites in the open, the recordings have to be scaled to the conditions inside the canopy (Gouttevin et al. 2015, Moeser et al. 2015. We used the inside-canopy modification of the meteorological conditions over the ground snow surface as realized in the hydroclimatological model AMUNDSEN (Strasser et al. 2011) and implemented the respective formulas in the new version of WaSiM. These are applied for the interpolated meteorological fields of solar and thermal radiation, air temperature, humidity and wind speed. The only parameter required is the effective leaf area index, LAI eff (-), e.g. from the literature, which is referred to herein as LAI for simplicity, i.e. including the area of the leaves or needles, the branches and the stem (Chen et al. 1997). Here, we present a condensed description of the canopy module as comprehensively described in Strasser et al. (2011).
Solar radiation reaching the ground surface Q sc↓ (W m -2 ) is calculated as a fraction of top-of-canopy incoming solar radiation Q s↓ transmitted through the trees depending on LAI (Hellström 2000): where 0.71 is a dimensionless extinction coefficient (Liston and Elder 2006). Longwave radiation reaching the ground Q lc↓ (W m -2 ) consists of a fraction of the top-of-canopy incoming longwave radiation Q l↓ , and a fraction of longwave radiation emitted by the trees: where σ is the Stefan-Boltzmann constant, T c (K) the canopy air temperature and F c (-) the canopy density, after Pomeroy et al. (2002), which was empirically derived from a comprehensive forest structure dataset: For the calculation of T c the dampening effect of the shading during the day and the emission of thermal radiation during the night are considered (Obled 1971): where T a (K) is the top-of-canopy air temperature, R c = 0.8, T mean (K) is the mean daily air temperature and δT is limited to the range -2 K ≤ δT ≤ + 2 K (Durot 1999): The increase of relative humidity RH c (%) inside a canopy due to sublimation and evaporation of melting snow is modified (Durot 1999): For melt conditions, RH c is set to saturation. Wind speed W c (m s -2 ) inside a canopy is calculated based on (Essery et al. 2003): where f i is the canopy flow index: with β (= 0.9) being a dimensionless scaling factor defined by Cionco (1978).

Canopy snow processes
For the computation of snow interception, sublimation as well as unloading by meltwater drip, and mass unload depending on effective LAI, we mostly followed the parameterization developed by Liston and Elder (2006) and Strasser et al. (2011). Again, we concentrated on the important formulas as implemented in the new WaSiM version. Snow that is intercepted on the leaves/needles, branches and stem of a tree can sublimate into the atmosphere, orif the air temperature is above 0°Cmeltwater drips and falls to the ground can occur. The unloading of snow by wind is not (yet) considered. Absorption of solar radiation SR abs (W) by snow in the canopy is given by: where the radius of a spherical ice particle r (m) is assumed to be 500 μm (Liston and Elder 2006); and Q s↓ (W m -2 ) is the top-of-canopy incoming solar radiation. For the intercepted snow particle albedo α, we assumed the same value as for simulated snow surface albedo in the open: where α min is the minimum albedo of (old) snow; α t-1 is the albedo in the previous time step; and k is a recession factor depending on air temperature (which determines snow surface temperature). The factor 1/24 is required to scale the result to the hourly progression of the computations. At every time step with a considerable snowfall (at least 0.5 mm h -1 ), the snow albedo is reset to its maximum value α max . For the calculation of the mass loss rate from the interception storage, the model requires the Reynolds, Nusselt and Sherwood numbers. The Reynolds number, Re, for 0.7 < Re < 10, is given by (Lee 1975): where v is the kinematic viscosity of air (1.3 × 10 -5 m 2 s -1 ). The Nusselt number, Nu, is given by: In contrast to the original version, we calculated the saturation vapour pressure over ice e s (Pa) after Alduchov and Eskridge (1997) Absolute humidity at saturation ρ v (kg m -3 ) was calculated using (Fleagle and Businger 1981): where the value of the gas constant for dry air, R d , is 287 J K -1 kg -1 . The diffusivity of water vapour in the atmosphere D v (m 2 s -1 ) was computed after Thorpe and Mason (1966): For the calculation of the mass loss rate (@m=@t) both temperature and humidity are assumed to be constant with height: where l s , the latent heat of sublimation, is 2.838 × 10 6 J kg -1 . In this computation, the Sherwood number, Sh, is set to Nu, and Ω is computed using: where λ t is the thermal conductivity of the atmosphere (0.024 J m -1 s -1 K -1 ); M w is the molecular weight of water (0.018 kg mole -1 ); and R is the universal gas constant (8.313 J mole -1 K -1 ). The sublimation loss rate coefficient ψ s (s -1 ) can now be computed using: where m sp (kg) is the particle mass (ψ i is the ice density = 916.7 kg m -3 ): The maximum snow interception storage capacity I max (mm) and the amount of solid precipitation in the current time step t, P (mm), were used to calculate the canopy-intercepted load I, given by : Rain is assumed to fall through to the rain interception model (to account for interception in the understorey) and finally to the ground snow cover. The maximum interception storage capacity I max is assumed to be 4.4 × LAI (Hedstrom and Pomeroy 1998), and the sublimation loss rate Q cs (mm) for the intercepted snow within the canopy is: Sublimation thus only occurs at the surface of the intercepted snow. This is accounted for with the nondimensional canopy exposure coefficient C e (Pomeroy and Schmidt 1993): The shape of the intercepted snow is described by the dimensionless coefficient k C = 0.01 (Liston and Elder 2006). For snow that is removed by snowmelt unload, it is assumed that initial melt at its surface is required. The rate of melt unload L m (kg m -2 ) is calculated using the temperature index melt model of Pellicciotti et al. (2005), accounting for shortwave radiation and albedo: where C t is the temperature factor (0.05 mm h −1 K −1 ), and C α the albedo factor (0.0094 m 2 mm W h −1 ). The unloading rate is limited to 5 kg m -2 d -1 K -1 , after Liston and Elder (2006). The unloaded mass is added to the ground snow cover beneath the trees. In order to distinguish liquid from solid unload, an empirical computation of the load of snow (solid fraction) falling to the ground is introduced (Liston and Elder 2006): In the extended WaSiM version, liquid outflow L m is redirected to the rainwater interception model in order to account for rainwater interception in the understory. The solid part L m,solid , however, is accumulated to the snow pack in the snow model. Validation of the snow interception model has been conducted by Montesi et al. (2004) in the Fraser Experimental Forest (39°53′N, 105°54′W) of the US Department of Agriculture (USDA). In the framework of the SnowMIP2 project, the set of parameterizations was compared to a variety of other approaches for snow-canopy interaction simulation (Rutter et al. 2009).
The new model extension includes three parameters that have been estimated a priori without any further consideration in the subsequent calibration procedure. A threshold value has to be provided for LAI and the roughness length z 0 , respectively. In our study, we set these values to LAI min ¼ 1.0 and z 0;min ¼ 0.1 m. Additionally, a scaling factor n LAI is mandatory for increasing the internal LAI values used for evapotranspiration calculation in WaSiM in order to consider the effective LAI for snow interception. The dense structure of needles and branches at the scale of a tree generally results in higher interception storage capacity for snow, exceeding typical values for rain Halldin 2001, Koivusalo andKokkonen 2002). LAI values provided in the lookup tables in the WaSiM settings (originally used for rain interception and evapotranspiration) are hence corrected by a scaling factor of n LAI ¼ 1.2 in order to account for the effects described above. This value leads to effective LAI values comparable with those in the literature for similar forests (Strasser 2008).

Model calibration and validation
The model was calibrated using a lexicographic calibration scheme introduced by Gelleszun et al. (2017). By means of this process-oriented approach, parameters are mapped to processes and an order of processes is suggested in which each is subject to a single calibration step. Here we performed two calibration steps in sequence: the first step includes parameters of the soil model (flow density and recession of saturated hydraulic conductivity with soil depth). This calibration step focuses on the calibration of interflow. The Nash-Sutcliffe model efficiency was the objective function used in this calibration step. The downhill simplex method (Nelder and Mead 1965) was applied to optimize two parameters simultaneously, according to Gelleszun et al. (2017). Likewise, in the second step, the parameters of the groundwater recession model were calibrated by minimizing the RMSE of the lowermost 35% of the flow duration curve (scaling parameter and recession constant). The subdivision into two single calibration steps with two different objective functions allows automatic calibration of the model with an acceptable number of iterations (100-150). This scheme was applied to both catchments in a similar way. All other model parameters were not calibrated, since for the physically-based energy balance snow model (Warscher et al. 2013) no calibration is required. Table 3 lists the parameters achieved by model calibration for the Brixenbach and Sieber catchments. Apart from the four parameters altered in the lexicographic calibration approach, only the precipitation correction was adjusted in a different way for each catchment, since first tests with the model suggested that the under-catch in the Sieber catchment is higher than that in the Brixenbach catchment. The calibration period is rather short due to data availability constraints; however, air temperature and precipitation recordings of these years show sufficient variability in terms of deviations from average conditions. For instance, in 2010 the mean annual temperature in the Brixenbach catchment was 4.2°C and mean annual precipitation amounted to 1500 mm. The respective values for 2011 are 6.0°C and 1450 mm, while 2012 was above average in terms of annual mean precipitation (5.2°C, 1950 mm). Thus, the calibration period reflects a broad range of climatic conditions.
For the Brixenbach catchment, the Nash-Sutcliffe model efficiency yields only 0.5 in the calibration period, and a percentage bias (PBIAS, see, e.g. Moriasi et al. 2007) value of about 10%, which is mostly related to mismatches in 2012 (this can be seen in the "Results and discussion" section, Fig. 6). The other parameters, including root mean square error (RMSE), RMSE to standard deviation ratio (RSR), and Pearson correlation coefficient (see Table 4), underline that the model tracks the observed time series reasonably well, given that only four parameters are adjusted in the calibration process. However, for the validation period, the model of the Brixenbach catchment shows better performance. The Nash-Sutcliffe model efficiency increases to 0.67, which can be considered as reasonable model skill. Similarly, the other performance measures indicate a better agreement of simulated with observed discharge valuesexcept for the RMSE, which increases slightly. However, the PBIAS, the RSR, and the Pearson correlation coefficient all show an increase in model skill in the validation period, which is in line with the higher Nash-Sutcliffe model efficiency value.
While for the Sieber catchment most skill measures reflect better performance compared to the Brixenbach catchment, the PBIAS and the RSME are higher in the Sieber catchment. Possible reasons might include incomplete knowledge of the actual operation of the Upper Harz Water Management System, which has been included in the model in a simple way using long-term averages, and the under-representativeness of rainfall stations in the higher elevation bands. However, the Nash-Sutcliffe model efficiency of around 0.8 suggests good model performance for both the calibration and the validation periods. The timing is also well captured by the model, which is obvious when considering the Pearson correlation coefficient of 0.9. The other measures show values similar to the respective values achieved for the Brixenbach catchment. Given that large portions of the catchments are steep slopes with shallow soils and with mostly unknown hydraulic characteristics of the bedrock and the lower model performance of hourly simulations (when Table 3. Parameter settings for both catchments. Parameters in italics were calibrated for each catchment using a lexicographic calibration approach. Parameters of the precipitation correction module were adjusted for each catchment. Each of the remaining parameters listed in the table is the same for both catchments.   Moriasi et al. 2007), the model set-ups can be considered as decently representing the water balance components required for an analysis of the effects of snow-canopy interaction on winter hydrology, as considered here.

Results and discussion
Snow cover dynamics Besides runoff, snow observations gathered in the Brixenbach catchment were also considered for assessing the model performance. This is especially relevant since in our study we wanted to focus on the influences of snow processes on streamflow. Figure 3 shows simulated sub-canopy SWE as well as intercepted snow load in the forest canopy in the Brixenbach catchment on 8 March 2016 (Fig. 3(a)) and 12 March 2016 ( Fig. 3(b)). The areas indicating a value of 0 (blue) in the canopy interception figures show the open areas of the study catchment. For the forested areas, the intercepted snow was between 16 and 20 mm on 8 March 2016. The values of sub-canopy SWE in the forest range between 14 and 147 mm, with increasing amounts at increasing altitude. That variation is also valid for the SWE values at the open site, but with higher amounts. While the SWE results for 12 March 2016 remain at comparable values, the snow interception load decreases and ranges from 0 to 13 mm on that day. Even though melt unload contributes to an increase in the SWE under the trees, the relative changes in canopy interception are higher. Figure 4 shows the observed and simulated SWE dynamics at the three study locations, Choralm, Doppelkehre and Talkaser Alm, in the Brixenbach catchment. The simulated values were extracted from the respective model raster cell. It becomes evident that the modified model is generally able to capture the dynamics of the snow cover at the observation locations well.
The accumulation and the ablation are well reproduced at the forest location at Choralm. There are differences in the peak (later and smaller peak in the simulations) and the duration of the snow cover is two days longer in the model simulation. The ablation period of the open location at Choralm was simulated very well. However, there was a systematic underestimation of the SWE amount by the model, which is especially obvious for the Talkaser Alm site. Possible reasons for this mismatch include underrepresentativeness of spatial precipitation variability and snow redistribution by wind at small spatial scale, both of which cannot be captured by the model.
The model performance for the SWE simulation of the open location at Doppelkehre is very good. Peak SWE was simulated with some delay and the duration of the snow cover simulated was only a few days too long, compared to the observation. The simulation of the snow cover at the adjacent forest site shows some limitations. The general dynamics and quantities were met fairly well, but the model was not able to simulate the complete melt and subsequent accumulation during the winter season. This limitation also became evident in SnowMIP studies at subalpine or maritime locations (Etchevers et al. 2004, Rutter et al. 2009) and represents a general pattern visible in the simulation results. Although the horizontal distances among the stations are very low (<2 km), the three sites show remarkable differences regarding snow dynamics during the winter season. This might be an indication for the fact thateven though the model is generally capable of representing the variability of processes at these scalesthe spatial variability is not covered by the meteorological input data. The usage of single precipitation stations (only one in the catchment and the others in the vicinity of the study catchment) and the 1-km-scale grids of all other meteorological variables might be a reason for these deficiencies.
Given that comparison between a point-scale observation and a grid cell of a model is always subject to mismatches in scale, and that the study catchment is very heterogeneous, the model results can be regarded as a good representation of the pattern and dynamics of the snow cover development. In principle, the differences between open and forested sites are captured surprisingly well. Figure 5 demonstrates how snow interception is computed by the model at the different sites of the Brixenbach catchment. Observed states are included in the plots as well in order to assess the accuracy of the model. In contrast to modelled time series, which represent SWE values intercepted by trees, the observed state extracted from the images is binary (snow on the trees or not). Comparison of the observed and modelled interception dynamics shows very good agreement. There are only a few events for which the model simulated snow in the forest canopy without evidence in the images and vice versa. Given that phase determination (liquid or solid) is always subject to uncertainties (Mair et al. 2016), the simulation results are very good, with just a few incorrectly modelled interception events.

Discharge simulations
Since discharge represents an integrative response of all hydrological processes in the catchment upstream, changes in representations of single processes can be smoothed out by the interplay of these processes at the catchment scale (Müller et al. 2017). Figure 6 shows the discharge simulated with the original and the extended WaSiM version compared to the observations at the outlet of both catchments. The simulations compare well with the observations. Some peaks in spring are lower in the extended model version when compared to the original model. This is in line with previous findings discussed for snow in the canopy. The snow under the canopy is decoupled from turbulence and radiation. This is one reason why the melt rates are typically smaller in the forest. This also explains the lower melt rates in the extended model version which, on average, shows better model performance in both forested catchments than the original version. However, in most periods, both model versions perform similarly.
In order to show the effect of the canopy model during periods subject to snowmelt, discharge in late winter and spring is plotted for both catchments in Figure 7. We selected the spring period 2015/16 during which the cameras were installed in the Brixenbach catchment. As this period was not available for the Sieber catchment, the rain-on-snow event in spring 2006 was selected, which corresponds to a return period of approximately 5 years. In general, the previous finding also holds true with respect to the simulation of snowmelt shown in Figure 7. With the extended model version, a lower snowmelt peak is computed, and thus a better match with the observed time series. In terms of performance measures, this effect is more relevant at the event scale. The Nash-Sutcliffe efficiency increases noticeably in both cases. When switching from the original to the extended version, the Nash-Sutcliffe efficiency increases from −0.68 to 0.71 (Brixenbach) and 0.21 to 0.87 (Sieber). This can be explained by reduced amounts of snow being present in the forested areas of the catchments prior to this melt event, and the lower melt rates due to the sheltering effect of the forest.
From Figure 8, which shows the corresponding basin-scale SWE computations for each catchment, it is obvious that the extended model computes not only lower accumulation of snow but also lower melt rates during spring. This is line with our findings discussed in connection with runoff during snowmelt in spring. Finally, Table 5 summarizes the comparison between both model versions in terms of the performance measures introduced for model calibration and validation. In contrast to the plots in Figure 7, the entire series are summarized in Table 5.
The comparison of the original (without canopy extension) and the extended model in terms of quantitative measures highlights that the model extension might help to better represent hydrological processes in forested subalpine regions with a seasonal snow cover. For both catchments the performance measures show improved values for the new model version that includes the canopy extension (see Table 5).

Conclusions and outlook
A new model extension is presented that makes the wellknown hydrological model WaSiM capable of representing snow interception and modified meteorological conditions in forests. Our work complements recent progress in enhancing the snow hydrology capabilities of the model (e.g. Daanen andNieber 2009, Warscher et al. 2013). Simulation results were compared to both observed discharge in two catchments and on-site automatic time-lapse photography, providing SWE and interception time series of snow at three different locations in one of the study areas. Even though the extended model version was calibrated and validated against discharge only without calibrating snow interception and snowmelt, the simulated intercepted snow and subcanopy snow accumulation match the on-site observations very well. These findings underline the importance of on-site observations for testing internal hydrological state variables. The extended model version also compares well with discharge observations, which has been shown for two catchments in different locations. Results suggest that the new snow-canopy interaction processes also prove to be helpful in simulating discharge.
The extended model introduces three additional parameters. Two of them, i.e. the thresholds for effective leaf area index (LAI) and roughness length z 0 , are only used to distinguish between forests and open sites. Hence, these values can be defined according to the land-use classification and are mostly insensitive to changes. Therefore, the LAI scaling factor n LAI is the only parameter needed for model tuning. With the LAI value chosen in this study, the model tracks interception, accumulation and melt in the forest very well. Possible mismatches in magnitude and timing of the seasonal cycle of snow interception and SWE accumulation are attributable to scaling issues introduced through comparing grid cells with point-scale observations and the limited spatial resolution of the meteorological forcing.
Even though the parameterizations utilized in this study are less complex than in similar models (Gouttevin et al. 2015), the good model reconstruction of observed on-site snow interception and SWE conditions in forests and on open sites suggests that the new WaSiM extension is suitable for regional-scale applications in hydrological modelling. Moreover, the introduction of only one additional model parameter for adjusting the LAI is regarded as a compromise between additional complexity and sophistication on the one hand, and a moderate increase in the degree of freedom for tuning the model on the other.
With the presented new model version of WaSiM it is now possible to accurately simulate the hydrologically very important and complex processes of snow accumulation and melt under forest canopies in a conceptually correct way. The significant value of the new forest parameterizations is given by the fact that they are part of a well-known, well-documented and worldwide applied hydrological model framework. The dissemination of the new model extension is foreseen in the next release of WaSiM and it will then be open source, hence freely available. were engaged in general support of the two projects and helped with the field work, image processing and figure preparation; Dr Stefan Pohl ( †) (University of Freiburg) participated in the development of the methodology of the time-lapse photography and helped with the analysis of the resulting digital images; Professor Dr Herbert Formayer and Dr Imran Nadeem (BOKU Vienna) provided the INCA data; Dr Jörg Schulla (Zürich) programmed the interface in WaSiM, and Dr Florian Hanzer (University of Innsbruck) helped to debug the new model extension; Alois Simon shared his local experience to interpret the soil maps, and Patricia Schrittwieser (both Federal State of Tyrol) helped in the parameterization of the tree species in the study area. Meteorological data were provided by the Hydrographic Service of Tyrol. Additional meteorological and hydrological data were kindly provided by the Harzwasserwerke GmbH, Hildesheim, Germany. Last but not least, the authors would also like to thank one anonymous referee, as well as Dr Michal Jenicek and the Associate Editor, Dr Guillaume Thirel, for their helpful comments and suggestions that significantly helped to improve the quality of the manuscript.

Disclosure statement
No potential conflict of interest was reported by the authors.