Estimating the snowfall limit in alpine and pre-alpine valleys: A local evaluation of operational approaches

The snowfall limit has important implications for different hazardous processes associated with prolonged or heavy precipitation such as flash floods, rain-on-snow events and freezing precipitation. To increase preparedness and to reduce risk in such situations, early warning systems are frequently used to monitor and predict precipitation events at different temporal and spatial scales. However, in alpine and pre-alpine valleys, the estimation of the snowfall limit remains rather challenging. In this study, we characterize uncertainties related to snowfall limit for different lead times based on local measurements of a vertically pointing micro rain radar (MRR) and a disdrometer in the Zulg valley, Switzerland. Regarding the monitoring, we show that the interpolation of surface temperatures tends to overestimate the altitude of the snowfall limit and can thus lead to highly uncertain estimates of liquid precipitation in the catchment. This bias is much smaller in the Integrated Nowcasting through Comprehensive Analysis (INCA) system, which integrates surface station and remotely sensed data as well as


A B S T R A C T
The snowfall limit has important implications for different hazardous processes associated with prolonged or heavy precipitation such as flash floods, rain-on-snow events and freezing precipitation. To increase preparedness and to reduce risk in such situations, early warning systems are frequently used to monitor and predict precipitation events at different temporal and spatial scales. However, in alpine and pre-alpine valleys, the estimation of the snowfall limit remains rather challenging. In this study, we characterize uncertainties related to snowfall limit for different lead times based on local measurements of a vertically pointing micro rain radar (MRR) and a disdrometer in the Zulg valley, Switzerland. Regarding the monitoring, we show that the interpolation of surface temperatures tends to overestimate the altitude of the snowfall limit and can thus lead to highly uncertain estimates of liquid precipitation in the catchment. This bias is much smaller in the Integrated Nowcasting through Comprehensive Analysis (INCA) system, which integrates surface station and remotely sensed data as well as outputs of a numerical weather prediction model. To reduce systematic error, we perform a bias correction based on local MRR measurements and thereby demonstrate the added value of such measurements for the estimation of liquid precipitation in the catchment. Regarding the nowcasting, we show that the INCA system provides good estimates up to 6 h ahead and is thus considered promising for operational hydrological applications. Finally, we explore the medium-range forecasting of precipitation type, especially with respect to rain-on-snow events. We show for a selected case study that the probability for a certain precipitation type in an ensemble-based forecast is more persistent than the respective type in the high-resolution forecast (HRES) of the European Centre for Medium Range Weather Forecasts Integrated Forecasting System (ECMWF IFS). In this case study, the ensemble-based forecast could be used to anticipate such an event up to 7-8 days ahead, whereas the use of the HRES is limited to a lead time of 4-5 days. For the different lead times investigated, we point out possibilities of considering uncertainties in snowfall limit and precipitation type estimates so as to increase preparedness to risk situations.

Introduction
The snowfall limit has important implications for different hazardous processes associated with prolonged or heavy precipitation. From a hydrological perspective, the snowfall limit in fact determines the percentage of liquid precipitation in a catchment and therefore exerts control on direct runoff (Tobin et al., 2012), as most precipitation at extra-tropical latitudes is formed from the solid phase (Fabry, 2015). In snow-covered areas, the presence of liquid precipitation is decisive for the occurrence of rain-on-snow events during which runoff can be considerably enhanced by snowmelt Morán-Tejeda et al., 2016;Würzer et al., 2016). Furthermore, in glaciated catchments, the aggregate state of precipitation at the surface will affect albedo and thus the production of glacial meltwater (Rohrer, 1989). For example, snow and glacial meltwater were relevant for flood generation in the Swiss Alps on 10 October 2011. During this event, the combination of meltwater and long lasting rainfall up to about 3000 m a.s.l. triggered debris flows and floods, which caused damage in the order of 94 million US$ or an equivalent of 71% of all precipitationinduced damage in that year (Andres et al., 2012). Characteristics of the snowfall limit are also relevant for the occurrence of freezing precipitation, which is not uncommon during winter months over continental Europe and northern America (Carrière et al., 2000;Cortinas Jr. et al., 2004;Rauber et al., 2001). Freezing rain, i.e. the occurrence of supercooled rain drops falling onto a sub-freezing surface (Forbes et al., 2014), is particularly dangerous because of its iceloading effects on power wires and potentially catastrophic cascading effects (Call, 2010;Chang et al., 2007). For example, after a prolonged episode of freezing rain in Slovenia in 2014, more than 300 power lines were broken and 25% of Slovenian residents were without electricity, heating and water (Forbes et al., 2014). Most frequently, freezing rain develops ahead of a surface warm front, where snow particles melt in an elevated warm layer and thereafter fall into a sub-freezing, nearsurface layer of very cold continental or Arctic air (Forbes et al., 2014;Rauber et al., 2001). Thereby, the depth of the melting layer determines whether the snow particles melt completely and turn into freezing rain, or whether they melt only partly and eventually refreeze to ice pellets until they reach the ground surface (Forbes et al., 2014;Reeves et al., 2014).
The snowfall limit can be estimated by operational approaches for different lead times, but its estimation remains challenging, especially in mountainous terrain (Schauwecker et al., 2016;Tobin et al., 2012). For hydrological applications, the snowfall limit is often estimated by an inter-and extrapolation of surface temperature measurements (Tobin et al., 2012). In those cases where temperature measurements are available within the catchment, extrapolations to higher and lower elevations are realized by assuming a constant lapse rate, which has been shown to be subject to considerable uncertainties (Schauwecker et al., 2016). If temperature measurements are, by contrast, not available within a catchment, a spatial interpolation approach is usually applied to more distant observations (Viviroli et al., 2009). In mountainous terrain, however, where the snowfall limit can be highly dynamic, such procedures will inevitably introduce additional uncertainty. In the case of steep valleys, the snowfall limit can be lowered locally by the cooling of the air during a single precipitation event (Unterstrasser and Zängl, 2006), a phenomenon which can hardly be predicted from distant observations. As an alternative to temperature inter-and extrapolation, the snowfall limit can also be inferred from radar observations, as they can detect falling hydrometeors directly and thus estimate their type from polarimetric properties of the radar signal (Chandrasekar et al., 2013;Ryzhkov and Zrnic, 1998). While such approaches have been shown to work well for stratiform precipitation over flat areas, they still suffer from limitations in the case of convective precipitation or if the radar signal is blocked by mountain topography at lower elevations (Besic et al., 2016). To predict the snowfall limit and precipitation type with a certain lead time, numerical weather prediction models are used; these are available for different temporal and spatial scales (Alfieri et al., 2012). However, direct outputs of such models may be of limited use for flood forecasting in steep terrain and in deep, narrow valleys, where model topography remains coarse and thus differs too much from actual conditions. Even worse so, in these areas, the height range with mixed precipitation can be broader than over flat areas under certain conditions (Haiden et al., 2011).
The accurate identification and forecast of freezing precipitation continues to be a substantial challenge for numerical weather prediction models (Forbes et al., 2014;Reeves et al., 2014). Recent improvements in the high resolution model (HRES) of the European Centre for Medium-Range Weather Forecast Integrated Forecasting System (ECMWF IFS) include a more representative timescale for the refreezing of rain drops, which is expected to improve the distinction between ice pellets and freezing rain (Forbes et al., 2014). Nonetheless, the generation of freezing drizzle drops from supercooled water in the absence of a warm melting layer is not yet represented in this system (Forbes et al., 2014). A detailed analysis of such an event and implications for the prediction of similar events can be found in Fernández-González et al. (2014,2015). Also, more experience needs to be gained with ensemble forecasts (ENS) of precipitation type, which are expected to increase the accuracy of mixed precipitation and in particular freezing rain forecasts (Forbes et al., 2014;Gascón et al., 2018).
Therefore, the aim of this study is to evaluate operational approaches estimating the snowfall limit and precipitation type for different lead times at the local scale. Thereby, we investigate the performance of temperature interpolation approaches as well as weather prediction from ECMWF IFS at different temporal and spatial scales in the Zulg valley, Switzerland. Resulting estimates are compared to measurements of a vertically pointing micro rain radar (MRR) and a disdrometer, which provide local information about these quantities at high temporal resolution. For the different lead times considered, the following questions are addressed: The paper is organized as follows: In Section 2, the study area is described and typical seasonal patterns of the snowfall limit are shown. In Section 3, available data and methods to estimate both the snowfall limit and related uncertainties are outlined. Results are presented in Section 4 and discussed in Section 5. Finally, conclusions with respect to operational monitoring and prediction of the snowfall limit are drawn in Section 6.

Study area
Operational approaches to estimate the snowfall limit are evaluated in the Zulg valley, Switzerland (Fig. 1). The Zulg catchment is 89 km 2 large and reaches from the confluence of the Zulg and the Aare river at 550 m a.s.l. to elevations up to 2060 m a.s.l. The bulk of the catchment (71% of the catchment area) is located between 800 and 1400 m a.s.l. After the confluence of three steep mountain torrents at 1030 m a.s.l., the Zulg descends to the outlet with a rather gentle mean slope gradient of 3%.
Estimation of the snowfall limit in this area is most relevant for hydrological applications during those periods of the year for which it is located within the elevational range of the catchment. The hypsography of the Zulg catchment is shown in Fig. 2 along with a climatology of the freezing level height at Payerne (see Fig. 3 for the location of Payerne), the latter being located typically 300-400 m above the snowfall limit (Fabry and Zawadzki, 1995). Between mid-April and mid-November, the median freezing level height is located above the highest mountain peaks of the Zulg valley. During this period, precipitation in the catchment can be assumed to fall predominantly in the liquid state. However, the freezing level height can strongly deviate from its climatological median. For example, a few days before the flood event in October 2011, it was located 1000-1500 m below its median, favouring conditions for a rain-on-snow event on 10 October 2011. Considering the 5th and the 95th percentile of observed freezing level heights, an accurate estimation of the snowfall limit in the Zulg valley is already important from the beginning of October and until the end of May. During this period, the snowfall limit can be located within the elevational range of the catchment and will thus determine whether most of the precipitation is falling in its solid or in the liquid state.

Data and methods
The operational approaches applied to estimate the snowfall limit in the Zulg valley are listed in Table 1 and the location of the different data sources considered is shown in Fig. 3. Information from the radio soundings in Payerne was used to determine large-scale seasonal patterns of the snowfall limit and the relevant period for its estimation in the Zulg valley. A temperature interpolation approach was further applied to monitor the snowfall limit by station measurements of both dry-bulb (T) and wet-bulb (T W ) temperature. Finally, products derived from numerical weather prediction models for both nowcasting and medium-range forecasting from ECMWF IFS are investigated. Whereas ECMWF HRES and ECMWF ENS are updated every 12 h, we only used updates every 24 h in this study. The local measurements of the MRR and the disdrometer and the data processing are then described before the metrics used to compare the operational approaches to these measurements.
Noteworthy, the snowfall limit typically refers to a fuzzy rather than a sharp boundary and its exact definition varies among different applications (Tobin et al., 2012). According to the American Meteorological Society (2017), we refer to the melting layer as the altitude throughout which solid precipitation melts as it descends during a precipitation event, and denote its upper and lower height limits as ZT and ZB, respectively. Furthermore, the freezing level height is defined as the lowest level in the free atmosphere for which the temperature is 0°C (American Meteorological Society, 2017).

Radio soundings
Information from radio soundings was used in this study to derive large-scale seasonal patterns of the snowfall limit and thus the relevant Median freezing level heights (1981-2010) as derived from radio soundings at Payerne (left) are compared to the hypsography of the Zulg valley (right). To illustrate the variability of the freezing level height during a specific year, deviations from the median are shown for 2011. Note that the snowfall limit typically lies 300-400 m below the freezing level (Fabry and Zawadzki, 1995). time period for its estimation in the Zulg valley. Two upper-air soundings per day (00:00 and 12:00 UTC) are available from the Payerne station of the Federal Office of Meteorology and Climatology (MeteoSwiss) since 1953, which allowed extracting a climatology of typical freezing level heights from 1981-2010 (Fig. 2). To extract the freezing level height from each radiosonde ascent, the vertical temperature profile was linearly interpolated between discrete measurements.

Temperature interpolation
In hydrological models, the snowfall limit is often determined from surface temperature measurements (Tobin et al., 2012). For example, in the European Flood Alert System (EFAS), averaged daily surface temperature measurements are interpolated using a temperature lapse rate of 0.65°C/100 m (Alfieri et al., 2014;Burek et al., 2013). Whereas most hydrological models consider only dry-bulb air temperature (Tobin et al., 2012), it has been shown that wet-bulb air temperature is a better index to discriminate between rain and snowfall, because air moisture influences the melting of snow (Froidurot et al., 2014;Rohrer, 1989). Therefore, we consider both dry-bulb and wet-bulb temperature measurements in this study. In Switzerland, such measurements are recorded every 10 min by the automatic monitoring network of MeteoSwiss (SwissMetNet). To estimate the temperature distribution in the Zulg valley from this dataset, the temperature interpolation approach described by Viviroli et al. (2009) was adopted. Thereby, both an elevation dependent regression (EDR) and a spatial interpolation of the resulting residuals (IDW) were applied to the measurements of 15 stations surrounding the study area, all located within distances < 35 km from the Zulg catchment and at elevations comprised between 471 and 3580 m a.s.l. (Fig. 3). Although we refer to this approach as temperature interpolation in the following, it should be noted here that its characteristics will be more similar to a temperature extrapolation outside this elevational range. To assess uncertainties inherent to this approach, analyses were performed independently for all 15 stations over the last three years (2015-2017) and results evaluated by crossvalidation. Since uncertainties are of interest only during rainfall events in this study, the validation period was restricted to those dates for which at least 30% of the selected stations recorded precipitation during the previous or following ± 3 h. As confirmed by the analysis, the temperature lapse rate during such periods is more stable and thus allows for a more accurate temperature interpolation than during periods without precipitation. For the same reason, inversion situations were identified with the help of the latest available radiosonde ascent at Payerne (Froidurot et al., 2014) and subsequently excluded from analysis.
From the interpolated temperature field, the aggregate state of precipitation (i.e. snow, mixed, rain) is usually determined by a threshold temperature or temperature range (Viviroli et al., 2009). For example, in EFAS, precipitation is assumed to be snow if the average daily temperature is below 1°C (Burek et al., 2013). Threshold temperatures for both dry-bulb and wet-bulb temperature were defined according to the approach of Froidurot et al. (2014). This approach estimates threshold temperatures from historical weather observations by estimating probability functions p rain and p snow in the form of Eqs.
(1) and (2), depending on temperature T and the two parameters α and β. Historical weather observations are used here from 15 synoptic stations from 1980-2017 ( Fig. 3), usually available three times a day since at least 1980, and at some stations even more frequently (e.g. Interlaken). These weather observations include information on the currently observed precipitation type (WMO4677), from which the aggregate state of precipitation was derived following Gjertsen Fig. 3. The location of different data sources considered to estimate the snowfall limit in the Zulg valley. Information from meteorological station is used for the application of temperature interpolation. The different data sources and their use are further described in Section 3.  Froidurot et al. (2014), only observations between −3 and 5°C were used to calibrate the models, as this has been proven to be the most relevant range to discriminate rainfall from snowfall.
When neglecting mixed precipitation, a threshold temperature discriminating rain from snow can be set at p rain = p snow = 50%. It should thereby be noted that the estimated probability will not correspond to the mixing ratio of rain and snow during a single event, but that it will rather represent the statistical likelihood of a temperature being associated with either an observed rainfall or snowfall event (Froidurot et al., 2014). In this study, the approach is extended by including mixed precipitation in order to define an upper and a lower limit of the melting layer. Thereby, p rain is calibrated by assuming all mixed precipitation to be snow, thus p rain = 50% is representing the lower boundary of the melting layer ZB. Accordingly, p snow is calibrated by assuming all mixed precipitation to be rain, thus p snow = 50% is representing the upper boundary of the melting layer ZT.
The resulting models and threshold temperatures to discriminate snowfall, rainfall, and mixed precipitation are shown in Fig. 4. Threshold temperatures associated with the upper and lower boundaries of the melting layer are estimated to be 0.8 and 1.4°C dry-bulb temperature and 0.3 and 0.8°C wet-bulb temperature, respectively. It can be seen qualitatively that the transition range between rain and snow is narrower for wet-bulb than for dry-bulb temperature and any consideration of the former will thus lead to a better discrimination of the two aggregation states, which is consistent with results from other studies (Froidurot et al., 2014).

The INCA system
Several parameters of the INCA system are provided by MeteoSwiss for the Zulg valley, namely the analysis and short-term predictions (i.e. 0-6 h) of the melting layer top (ZT) and the melting layer half width, from which the melting layer bottom (ZB) can be calculated. This system has a horizontal resolution of 1 km and was developed for applications in mountainous terrain, where direct outputs of numerical weather prediction models may be of limited use due to their limited horizontal resolution. The analysis of the INCA system combines surface station data with remotely sensed information in such a way that the observations at the station locations are reproduced, whereas the remotely sensed data provide the spatial structure for the interpolation. In the forecast, measurements are further combined with numerical weather prediction model outputs. To this end, the Swiss version of INCA uses the Consortium for Small-Scale Modeling (COSMO) fields (Haiden et al., 2011). The distinction between rain and snow is based on the vertical profile of the wet-bulb temperature at each grid point, derived from 3D temperature and humidity fields. The melting layer top and half width are assumed to correspond to wet-bulb temperatures of 0.2 and 1.3°C, respectively. The INCA system has been validated based on 174 stations over the eastern Alpine region by Haiden et al. (2011). The authors conclude that temperature is analysed with an average accuracy of 1-1.5°C, with the smallest errors occurring in lowland areas and at exposed mountain stations. In Alpine valleys, by contrast, the mean absolute error (MAE) is larger (1.5-2.5°C), mostly due to the lack of precise information about inversion heights (Haiden et al., 2011). The operational INCA verification in Switzerland shows a similar evolution with forecast times as obtained by Haiden et al. (2011): For ground stations, an average MAE of 0.4-0.6°C at analysis time increasing with time up to 1.5-2°C at 6 h (Matteo Buzzi, personal communication, December 22, 2017).

ECMWF global models
In this study we are using the precipitation type forecast from HRES and the probability of precipitation type ENS product (Gascón et al., 2018) from ECMWF IFS. The HRES precipitation type describes the type of precipitation at the surface. A precipitation type is assigned wherever there is a non-zero value of precipitation in the model output field. This variable is available in IFS since 2014 and provides 6 different precipitation types with a horizontal resolution of 9 km: rain, snow, wet snow, mixture of rain and snow, ice pellets, and freezing rain. The three last types will be considered as mixed precipitation in the development of this study. The melting and refreezing parametrizations in the IFS model are formulated in terms of the two prognostic variables for precipitation: rain and snow. Precipitation type is diagnosed from the ratio of rain and snow at the surface and the profile of precipitation and temperature above. Whether the precipitation refreezes or not in the cold air below also affects the temperature profile, resulting in colder temperatures in the layer if the particles remain as supercooled water. The latent heat of fusion is instead transferred to the surface with the rain freezing on impact and a relative warming of the surface and nearsurface temperature. The complete description of the algorithms that define the different precipitation types can be found in Forbes et al. The probability of precipitation type is a variable from ECMWF IFS based on ENS forecast (approximately 18 km horizontal resolution) and it calculates the probabilities of each precipitation type at each location. Also, this product provides much more detail regarding the probabilities for different precipitation types and includes information pertaining to the instantaneous total precipitation rate, which can be key for determining the severity of impacts, or, for example, potential freezing rain or snowfall events (Gascón et al., 2018). Three different categories of precipitation type are defined depending on the precipitation rate, from a minimum precipitation rate for each precipitation type to 0.2 mm h −1 , from 0.2 to 1 mm h −1 and greater than 1 mm h −1 . A minimum precipitation rate threshold is applied for each precipitation type. The main objective was to enforce a zero frequency bias comparing with a 4-month verification training period that uses present weather reports from surface synoptic observations. Minimum precipitation rates of 0.05 mm h −1 were assigned as the most suitable limit for snow and freezing rain forecast, 0.1 mm h −1 for rain mixed with snow and 0.12 mm h −1 for rain. In this study, we applied the probabilities of precipitation type based on ENS to reduce misses and false alarms.
The last point to consider in this study was the difference between the height used in the model and the real height in the field (Gjertsen and Ødegaard, 2005). To deal with this issue, we considered the grid point being nearest to the study location and showing the most similar height with respect to the location of the instrumentation. For both the HRES and the ENS forecasts, this point is located at a distance of 25 km.

Local measurements
Since early 2017, a vertically pointing MRR and a disdrometer are installed in the Zulg catchment at 1060 m a.s.l. (Fig. 3). These devices provide information about the snowfall limit and precipitation type at the ground surface and with high temporal resolution (1 min). An example of the information gathered from these two instruments is shown for late April 2017 in Fig. 5. Here, the snowfall limit was decreasing in the atmosphere over the study area, causing a transition from rain to snowfall at the ground surface. The processing of the raw data from both instruments is described in the following.
The MRR used here is a 24-GHz (K-band) frequency-modulated continuous wave (FMCW) Doppler radar and was installed on 21 April 2017. It has a small transmission power (50 mW), such that a common antenna is used to transmit and receive power and that no beam overlapping occurs (Peters et al., 2005). The MRR captures the Doppler spectra at each gate from which the raindrop size distribution (DSD) is retrieved under certain assumptions. Also, the two-way attenuation from the scattering volume is estimated from the retrieved DSD to adjust the 24-GHz returned power into a Rayleigh scattering equivalent reflectivity factor (Tokay et al., 2009). In this study, a resolution of 100 m for each of the 30 gates was chosen, restricting the measurements to 3000 m above the ground surface (i.e. 1060-4060 m a.s.l.). From the measured Doppler spectra, ZT and ZB were determined according to the bright band detection algorithm described by Pfaff et al. (2015). Thereby, ZT and ZB are determined by the largest positive and negative reflectivity gradients, whereas the vertical profile of Doppler velocity is considered as an additional constraint. Furthermore, several filters are applied to the original measurements. First, only periods for which rain was detected by the disdrometer at the ground surface were considered, which is a basic condition for the melting layer to be observable by the MRR. Second, several plausibility checks were conducted with respect to the melting layer width (Schauwecker et al., 2017). More precisely, values were only considered plausible if the melting layer width could be defined (i.e. present upper and lower limits) and was within a range of 0-500 m. Also, profiles with a mean reflectivity below 10 dBZ were considered to be without bright band (Pfaff et al., 2015). Third, outliers beyond one standard deviation of ± 60 min were excluded, and the resulting values were smoothed over this time period. An example of raw reflectivity profiles as well as the melting layer obtained by the described algorithm is presented in Fig. 5.
The disdrometer measures the attenuation of a laser beam (785 nm) by precipitation particles with an optical sensor and was installed on 3 February 2017. From the amplitude and duration of this attenuation, the diameter and fall velocity of precipitation particles are derived and translated into precipitation type by considering the statistical relation between these quantities (Gunn and Kinzer, 1949). Given the estimated precipitation type, the aggregate state of precipitation (i.e. snow, mixed, rain) was derived for this study by using the same reclassification as for the synoptic weather observations described in Section 3.2.

Metrics
For the monitoring and the nowcasting range, we quantitatively compare estimates of operational approaches (x ) to MRR and disdrometer observations (x). In a first step, we characterize the bias of these approaches to assess potential systematic errors in estimates. For the continuous variables of ZT and ZB, which are compared to MRR observations, we thereby calculate the bias (B) and the correlation coefficient (Corr) according to Eqs. (3) and (4). For the multi-categorical variable of precipitation type (i.e. snow, mixed, rain), which is compared to disdrometer observations, we calculate the frequency bias B k for each category k, i.e. the ratio of the number of forecasts of 1 occurrence to the number of actual occurrences (Jolliffe and Stephenson, 2011).
In a second step, we characterize uncertainties after a bias correction based on MRR measurements. Thereby, we adjust the estimates of ZT and ZB to these measurements by linear regression for all considered approaches, and characterize the scatter of the remaining estimates. For the continuous variables of ZT and ZB, we thereby calculate the MAE according to Eq. (5). For the variable of precipitation type (i.e. snow, mixed, rain), we calculate the Gerrity equitable score (GS), as this score utilizes off-diagonal information in the contingency table and penalizes larger errors more (Jolliffe and Stephenson, 2011), e.g. in the case of observed rainfall, the forecast of snowfall will be more penalized than the forecast of mixed precipitation.
The medium-range forecast of precipitation type is investigated on the basis of a case study. We thereby focus on the comparison of the HRES and the probability of precipitation type ENS product from ECMWF IFS. To assess the temporal persistence of these forecasts according to Thielen et al. (2009), we visualise the observed as well as the predicted aggregate state of precipitation for different lead times in a colour-coded box representation.

Monitoring
The monitoring of the snowfall limit by surface temperature interpolation is compared to the analysis of the INCA system. Thereby, uncertainties related to surface temperature interpolation are assessed on a regional scale, before estimates of the snowfall limit and precipitation type are compared to local observations in the Zulg valley.
Residuals of the applied temperature interpolation approaches are shown in Fig. 6, as obtained by cross-validation. The interpolation of wet-bulb temperatures generally provides better results (MAE: 0.54°C) than the interpolation of dry-bulb temperatures (MAE: 0.78°C). The relevance of these uncertainties for the estimation of the snowfall limit depends on the lapse rate in the atmosphere and is less critical in situations with larger lapse rates. Given mean modelled lapse rates of − 0.58°C/100 m (dry-bulb) and − 0.55°C/100 m (wet-bulb) in the validation period, we would expect a MAE in snowfall limit estimates of 135 m (dry-bulb) and 100 m (wet-bulb), respectively. For the INCA system, assuming typical MAEs of 0.4-0.6°C for dry-bulb temperature at analysis time (Section 3.3), we would expect a MAE in snowfall limit estimates of about 70-105 m.
When comparing estimates of the snowfall limits to local MRR measurements, uncertainties are generally higher than in temperature estimates and a bias with regard to ZT and ZB is apparent in all the considered approaches. As can be seen in Fig. 7 (left) and Table 2, the applied temperature interpolation approaches seem to systematically overestimate the altitude of the snowfall limit, especially its lower limit ZB. When compared to disdrometer measurements, these approaches consequently overestimate the occurrence of mixed precipitation at the expense of snowfall at the ground surface. Also, differences between the interpolation of dry-bulb and wet-bulb temperatures become very small. More accurate estimates are provided by the INCA system. Here, a smaller positive bias is observed for ZT, whereas a slight negative bias is obtained for ZB. Thus, the occurrence of mixed precipitation is overestimated at the expense of rainfall at the ground surface, which is confirmed by the comparison to disdrometer observations. Noteworthy, the overestimation of mixed precipitation can generally be expected from deterministic approaches based on temperature thresholds: As shown in Fig. 4, a certain probability of pure rain-or snowfall between the assumed threshold temperatures will always exist, a fact which is neglected by all approaches. Because the INCA system assumes mixed precipitation to occur within the largest temperature range, the overestimation is also largest in this system (Table 2).
Furthermore, Fig. 7 shows that the bias seems to depend on both elevation (left) and the hour of the day (right). In particular, the interpolation of surface temperatures tends to strongly overestimate the altitude of the snowfall limit around noon. This is, the mean bias of ZT and ZB between 08:00-12:00 UTC rises up to about 300 and 500 m, respectively. This particular overestimation results from the exaggerated daily cycle of the snowfall limit inherent to such approaches and is much less apparent in the INCA system. The latter takes the three-dimensional properties of the atmosphere more accurately into account, therefore resulting in a much smoother daily cycle of estimated snowfall limits.
After applying an elevation dependent bias correction with respect to MRR measurements, uncertainties of the considered approaches are reduced. Remaining uncertainties are characterized in Fig. 8 and Table 3. Remaining MAEs of ZT and ZB are in good agreement with the expectations stated above, i.e. they amount to 126 and 125 m for drybulb temperature interpolation, 113 and 111 m for wet-bulb temperature interpolation, and 88 and 96 m for the INCA system, respectively. To assess the relevance of the described bias correction for the estimation of liquid precipitation in the catchment, we have investigated a specific case that occurred between the April 25, 2017 12:00 and April 28, 2017 12:00 UTC, when the snowfall limit was slowly decreasing during a precipitation event. In a first step, the precipitation type (i.e. snow, mixed, rain) is estimated every 10 min for the entire catchment area based on a temperature interpolation (T/T W ) and the INCA system. In a second step, the elevation dependent bias corrections based on the MRR measurements was applied to each of these approaches, before the resulting precipitation type was estimated once again. For each time step, we then calculated the catchment area for which these two estimates differ, so as to delimit those areas for which precipitation type estimates are uncertain. On average, temperature interpolation resulted in an uncertain prediction of precipitation type in 20% of the catchment area (using both T and T W ), whereas uncertainty was reduced to 12% when using the INCA system. Whereas these values reflect mean uncertainties during the entire event, they can be much larger at particular points in time. Fig. 9 shows the estimated precipitation type by interpolation of dry-bulb temperature with and without an application of the bias correction at the moment of maximal sensitivity, i.e. when precipitation type estimates differ in 45% of the catchment.

Nowcasting
To assess the potential of the INCA system for the nowcasting range, we compare predictions of this system for lead times of up to 6 h with local MRR and disdrometer measurements. The most relevant metrics are presented in Table 3 and complemented with their values for the monitoring range for reference. All metrics are calculated after an elevation dependent bias correction based on the MRR measurements. Table 3 and Fig. 10 demonstrate that the predictions of ZT and ZB provided by the INCA system with a lead time of up to 6 h perform better than the real-time estimates derived from surface observations. When comparing predictions of precipitation type to disdrometer observations, the skill of the INCA system is comparable to real-time estimates based on temperature interpolation up to a lead time of 6 h. When taking a closer look at different series of lead times, it can be seen that the skill of the prediction seems to decrease for the first 3-4 h and,  Fig. 7. Modelled upper and lower melting layer height limits (ZT and ZB) as derived from interpolation of dry-bulb (T) and wet-bulb (T W ) surface temperatures as well as the INCA system are compared here to MRR measurements (21.4.2017-31.12.2017) (left). The comparison is made before an elevation dependent bias correction. Residuals are shown for different hours of the day (right).

Table 2
Characterization of the bias in temperature interpolation approaches (T and T W ) and the INCA nowcasting system (lead times 0-6 h) compared to local observations of a micro rain radar (21.4.2017-31.12.2017) and a disdrometer (3.2.2017-31.12.2017 inversely, to increase again afterwards in the metrics MAE ZT and GS. In the MAE ZB , however, this pattern is not apparent.

Medium-range forecasting
The medium-range forecasting from ECMWF probability of precipitation type product and HRES precipitation type variable are explored for the snowfall limit on the basis of a selected case study between March 7-10, 2017. The precipitation type observed by the disdrometer at 1060 m a.s.l. during this period is given in Fig. 11 (bottom line). On March 7, 2017, and on the days before, the snowfall limit was located between 500 and 1000 m a.s.l., i.e. around the lower elevational range of the catchment. Precipitation during this period can be assumed to be snowfall within most of the catchment. Between March 8 and 9, 2017, the snowfall limit sharply increased from below 500 to above 2000 m a.s.l. This increase was accompanied by 50 mm of rainfall measured by a pluviometer at 1060 m a.s.l. Peak discharge of the Zulg river reached 35 m 3 s −1 , corresponding to the second highest discharge within the observation period (March 2-September 30, 2017). It can thus be assumed that the combination of the melting snow and the rainfall associated with the rise of the snowfall limit was important for runoff generation in this situation.
Predictions of precipitation type by the HRES and the probability of precipitation type from the ENS of the ECMWF are shown in Fig. 11 for the available lead times. For HRES, each box is coloured according to the predicted precipitation type on the ground. For the ensemble-based forecast, if precipitation is forecasted, each box is coloured proportionally to the probabilities of the predicted precipitation types. The rain event of March 9, 2017, is predicted by the HRES with quite good temporal persistence 4-5 days ahead of the event. However, for larger lead times, the forecast changes from one day to the next between precipitation and no precipitation (e.g. March 3-4, 2017) or even between different aggregate states of precipitation (e.g. Feb 28-March 1, 2017). The ensemble-based product, by contrast, shows a more persistent pattern at these lead times. Already prior to March 1, the snowline was predicted to be around this critical elevation and that a possibility for a rain-on-snow event could emerge at the study site. INCA T W T Fig. 8. Modelled upper and lower melting layer height limits (ZT and ZB) as derived from interpolation of dry-bulb (T) and wet-bulb (T W ) surface temperatures as well as the INCA system are compared here to MRR measurements (21.4.2017MRR measurements (21.4. -31.12.2017. The comparison is made after an elevation dependent bias correction. Residuals are shown for different hours of the day (right).

Table 3
Characterization of the scatter in temperature interpolation approaches (T and T W ) and the INCA nowcasting system (lead times 0-6 h) compared to local observations of a micro rain radar (21.4.2017-31.12.2017) and a disdrometer (3.2.2017-31.12.2017). All metrics are calculated after an elevation dependent bias correction. persistent in the forecast. This is also the case on March 4, where the HRES predicts no precipitation at all. Although the rain-on-snow event is associated with quite low probabilities in the ensemble-based forecast, the possibility of this event is visible more in advance and with more persistence than in the HRES.

Discussion
In the following, implications of the results for the estimation of the snowfall limit and precipitation type are discussed for the different lead times considered.
The monitoring of the snowfall limit in hydrological modelling is often based on surface temperature measurements of meteorological stations (Tobin et al., 2012). Several studies have shown that wet-bulb temperature is a better index to discriminate between rain and snowfall than dry-bulb temperature, because air moisture influences the melting of snow (Froidurot et al., 2014;Rohrer, 1989). However, a regionalization of air moisture has only rarely been investigated in the past (Froidurot et al., 2014). In this context, we demonstrate that the interpolation of wet-bulb temperature from distant observations (MAE: 0.54°C) is favourable to the interpolation of dry-bulb temperatures (MAE: 0.78°C). Furthermore, temperature interpolation approaches perform better during precipitation events, leading to lower interpolation errors as compared to studies where this restriction is not being applied (Jabot et al., 2012). However, when comparing estimates of the snowfall limit to local measurements of a MRR and a disdrometer, the interpolation approach seems to systematically overestimate the altitude of the snowfall limit and the occurrence of mixed precipitation instead of snowfall at the surface, respectively. This overestimation can be due to different reasons: (i) air temperatures measured 2 m above the ground surface deviate from air temperatures in the free atmosphere at the same elevation. In particular, the daily temperature cycle is more pronounced close to the surface than at higher elevations in the free atmosphere; (ii) in case of local precipitation and cloud cover, or in case of a local cold surface layer, temperature interpolation from distant observations will cause overestimations; (iii) during a single precipitation event, cooling by evaporation and fusion will cause a local decrease of the snowfall limit (Minder and Kingsmill, 2013 Table 3 describing the performance of the INCA system for lead times up to 6 h ahead. Estimates of the snowfall limit are compared to MRR measurements (i.e. MAE), while estimates of precipitation type are compared to disdrometer measurements (i.e. GS). We also indicate the performance of corresponding real-time temperature interpolation as a reference. All metrics are calculated after the elevation dependent bias correction. Unterstrasser and Zängl, 2006), which cannot be captured by distant observations.
In this study, we have evidence for the first of these reasons. More specifically, we show that temperature interpolation is not found to be biased if validated against surface temperature measurements, but seems to be biased if compared to direct measurements of the snowfall limit. In addition, this bias seems to increase slightly with increasing elevation. Our results are thus in agreement with findings from Schauwecker et al. (2017), who compared measurements of the freezing level height as obtained with a MRR to a simple temperature extrapolation approach from a single meteorological station at ground surface. Here, we complement these findings by showing that the overestimation is largest around noon, i.e. when warming of the surface layer is most pronounced. Finally, we also demonstrate that this bias is partly reduced by the INCA system, because it takes account of the three-dimensional properties of the atmosphere and because daily temperature variations are dampened at higher elevations in the free atmosphere. We thus also confirm results of Tobin et al. (2012), who stated that the use of the COSMO model will result in better estimates of the snow covered area and discharge than conventional interpolation of ground measurements. Despite these improvements, MAEs of the estimated melting layer height limits in INCA remain at values of 131 and 118 m for ZT and ZB, respectively. Such uncertainties can be especially relevant whenever the snowfall limit is located at critical elevations with respect to catchment hypsography, which is typically the case in the Zulg catchment from early October to late May. To reduce these uncertainties further, an elevation dependent bias correction was applied here based on local MRR observations. We show that the consideration of the hour of the day can reduce the particularly large overestimations around noon. In addition, and in the case of critical situations, a consideration of several scenarios would allow accounting for the variability of the melting layer width and related changes in the percentage of liquid precipitation in the catchment.
In the nowcasting range, we investigated predictions of the snowfall limit by the INCA system with lead times up to 6 h. We conclude that these predictions perform even better than commonly applied real-time estimates based on dry-bulb temperature interpolation, and that they should thus be regarded as a promising tool for local hydrological applications. Depending on the chosen skill score, we further find a slight increase in performance beyond lead times of around 3-4 h (i.e. for MAE ZT and GS). This behaviour is contradicting findings of Haiden et al. (2011), who reported a monotonous decrease in temperature forecast skill of the INCA system with increasing lead times, based on a validation using 174 stations over the eastern European region. However, in contrast to their study, we assessed estimates based on wet-bulb temperature profiles in the free atmosphere and compared them to direct observations of the snowfall limit provided by a MRR.
In terms of medium-range forecasting of the snowfall limit, we compared an ensemble-based product with the HRES of ECMWF to detect critical changes in the aggregate state of precipitation with lead times of several days. On the basis of a case study, we illustrate that an ensemble-based forecast has the potential to anticipate rain-on-snow events with lead times of up to 7-8 days. Although precipitation cannot be skillfully forecasted more than 2-3 days in advance (Thielen et Fig. 11. Medium-range forecasts of precipitation type during an observed rain-on-snow event in the Zulg valley. Disdrometer observations (OBS) are depicted in the lowermost line of each panel. Forecasts are depicted for different lead times and two products of the ECMWF: (I) The high resolution forecast (top), which is available for up to 10 days ahead and which predicts precipitation type deterministically, and (II) an ensemble-based forecast (bottom), which is available for up to 15 days ahead and which provides probabilities for different precipitation types. Here, each box is coloured proportionally to the probabilities of forecasted precipitation types.
2009), the forecast of temperature fields in relation with associated precipitation type is regarded here as a promising avenue for flood early warning. We thus agree with findings of EFAS, although established primarily for larger river basins (Demeritt and Nobert, 2014;Thielen et al., 2009), who concluded that ensemble-based systems can be used effectively to increase preparedness and internal vigilance to potentially critical situations at these lead times. This is even more so true as false alarms associated with a small probability of event occurrence are less critical to users than missed events ). Moreover, we show that the ensemble-based approach predicts the observed change in precipitation type more persistently than the HRES of ECMWF. Temporal persistence is an important criterion for the reliability of a forecast ) and will facilitate decision making. Therefore, and although the interpretation of ensemble-based forecasts for operational applications remains challenging (Demeritt and Nobert, 2014), our measurements support the use of such forecasts to issue pre-alerts with lead times of several days ahead.
Finally, it should be noted that for a given estimate of the snowfall limit and precipitation type at the ground surface, any prediction of whether precipitation falling in the form of snow or rain is being stored in an (antecedent) snow cover, evaporating, or contributing to runoff is not straightforward. Lundquist et al. (2008) investigated some cases for the Sierra Nevada (California) and offer possible solutions for this estimation.

Conclusions
In this study we have characterized uncertainties in estimating the snowfall limit for different lead times based on local measurements in the Zulg valley, Switzerland. When monitoring the snowfall limit, commonly applied temperature interpolation approaches seem to overestimate the altitude of the snowfall limit and underestimate the frequency of snowfall at ground surface. Whenever the snowfall limit is located at critical elevations within the catchment, such uncertainties may be relevant for hydrological modelling and operational flash-flood or debris-flow early warning systems. For local hydrological applications, we therefore suggest an elevation dependent bias correction based on local MRR measurements. Also, we provide evidence that this bias can be further reduced by taking into account the hour of the day to avoid particular overestimations of the snowfall limit around noon. In terms of nowcasting, the INCA system seems to perform better than real-time dry-bulb temperature interpolations for lead times of up to 6 h and is thus regarded as a promising tool for local operational applications. For the medium-range forecast, we conclude that an ensemblebased product of ECMWF has the potential to anticipate critical changes in precipitation type in the study area with lead times of several days. In our case, this product provided an earlier and more persistent signal of such changes than the HRES of ECMWF. However, further evaluation of ensemble-based forecasts of precipitation type is required to confirm these results. In this study, local measurements of the snowfall limit and precipitation type proved valuable for the evaluation of operational approaches at different temporal and spatial scales and the assessment of their suitability for local applications.