Analysis and Predictability of the Hydrological Response of Mountain Catchments to Heavy Rain on Snow Events: A Case Study in the Spanish Pyrenees

From 18 to 19 June 2013, the Ésera river in the Pyrenees, Northern Spain, caused widespread damage due to flooding as a result of torrential rains and sustained snowmelt. We estimate the contribution of snow melt to total discharge applying a snow energy balance to the catchment. Precipitation is derived from sparse local measurements and the WRF-ARW model over three nested domains, down to a grid cell size of 2 km. Temperature profiles, precipitation and precipitation gradient are well simulated, although with a possible displacement regarding the observations. Snowpack melting was correctly reproduced and verified in three instrumented sites, and according to satellite images. We found that the hydrological simulations agree well with measured discharge. Snowmelt represented 33% of total runoff during the main flood event and 23% at peak flow. The snow energy balance model indicates that most of the energy for snow melt during the day of maximum precipitation came from turbulent fluxes. This approach forecast correctly peak flow and discharge during normal conditions at least 24 h in advance and could give an early warning of the extreme event 2.5 days before.


Introduction
From June 18 to June 19, 2013, the Ésera river in the Pyrenees, Northern Spain, caused widespread damage due to flooding. Damage to public properties, such as roads or bridges, and to private property was estimated in excess of seven million euros. Over 300 people had to be evacuated from their homes, but fortunately there were no fatalities. The increased flow came as a result of torrential rains at a time when the catchment presented an anomalously extensive snow cover above 2000 m a.s.l. [1], and the 0 • C isotherm was above 4000 m a.s.l.
The role of rain on snow events (ROS) is interesting from a scientific and applied point of view, since it is widely recognised that they dominate much of the flood generation in mountainous and boreal regions [2,3]. Most of the largest floods in British Columbia, Washington, Oregon and California have been associated with ROS [4][5][6]. In Germany, ROS have much larger hydrological influence than rain alone in catchments above 400 m a.s.l. ROS have been identified as a primary cause of changes in channel morphology due to erosion [7][8][9][10][11]. For example, at a site in the Oregon Cascades, 85% of all landslides which could be accurately dated were associated with snowmelt during rainfall [7]. Sandersen et al. [12] reported that triggering of debris flows in Norway is caused by the combination of rainfall and snow melt. ROS events are also considered a major cause in the release of some type of avalanches [2].
The importance of ROS events is highlighted above, however there is uncertainty about the physical processes that control the runoff generation during their occurrence and their sensitivity to temperature and other meteorological variables. McCabe et al. [6] made a comprehensive study of several ROS events in the Pacific Northwest of the United States (PNW) concluding that the severity of a ROS events does not only depend on the magnitude of the precipitation but also on the elevation of the freezing level and the extent and characteristics of the existing snowpack. Marks et al. [13] studied a heavy ROS event (winter 1996) in the Pacific North West of the United States, where they found that 60-90% of the energy for snowmelt came from turbulent heat exchange. Condensation during the event was a significant contributor to the flood. Van Heeswijk et al. [14] concluded that wind speed together with vapour pressure and temperature gradients are the most important climate variables that control snow melt. However, in a later study Mazurkiewicz et al. [15] found that radiation was the largest contributor to melt during ROS in the Pacific North West. Other studies show the importance of the low permeability of the snowpack when it is saturated, leading to a quick generation of excess runoff. This indicates that heavy rain water moves several times faster than the natural snowmelt [2,11,16]. During ROS events the potential for flooding is increased if the soil is frozen [17,18].
Quantifying the actual contribution of snow melt and rain to runoff in a particular ROS event is challenging, as well as it is the extrapolation of results to other events or geographical areas. Such complexity makes it difficult to anticipate the hydrological response of a catchment, even when detailed meteorological forecasts are available. In the last decades the spatial resolution and the accuracy of atmospheric mesoscale models have improved rapidly. These models are now able to produce reliable meteorological fields over complex terrain when forced with realistic initial and boundary conditions [19][20][21]. Such fields can be used as inputs for hydrological models enabling the forecasting and analysis of extreme weather events, such as ROS, and their hydrological consequences. However, a small deviation in determining initial conditions, or biases in the input data for snow models may lead to considerable errors in the timing, magnitude and spatial distribution of the simulated flood [22,23].
In this study we investigate the skill of a high resolution weather model -the weather research and forecasting model, advanced research WRF core, WRF-ARW [24] -to simulate and forecast the occurrence of this extreme precipitation event. The output of the WRF combined with land observations and remote sensing data are used to feed a physically-based snow energy balance model. This approach is useful to improve our understanding of ROS in several aspects: 1. To assess the suitability of snow models coupled with atmospheric models to simulate ROS events 2. To estimate the spatial distribution of discharge in several sub-catchments 3. To estimate the relative forcing of meteorological variables to snow melt 4. To estimate the relative contribution of snowmelt and rainfall to the final river discharge 5. To estimate how much time in advance a weather model can give an early warning of dangerous floods All these questions are key to understanding the nature and predictability of ROS events in the Pyrenees and other similar mountain regions.

Site of Study
The Ésera river is one of the most important tributaries of the Ebro river in North Eastern Spain. Most of its runoff is generated in the headwater area located in the central Spanish Pyrenees. The flood event described in this study affected mostly the upper portion of the catchment, therefore, we restrict ourselves to the area above the gauging station of Eriste, at 1050 m a.s.l. (Figure 1). It was near this location that the most destructive effects of the flood were experienced. The catchment area is 323 km 2 and contains five small sub-catchments: Eriste, Estós, Maladeta, Ramascaró and Vallibierna. They drain one of the highest ranges in the Pyrenees, with many peaks exceeding 3000 m a.s.l., including the highest summit in the whole mountain range: Aneto (3404 m a.s.l.). In this catchment 65.6% of the area is located above 2000 m a.s.l., and 26.4% over 2500 m a.s.l. [25]. The geology exhibits a rather complex structure with granites, slates, schist and limestones heavily fractured and folded [26]. There is a mismatch in the extent of the topographic and the hydrological catchments, since several sinkholes divert part of the flow to the Garona basin (French Pyrenees). These diverted flow could represent 12-25% of the annual runoff of the catchment [25,27]. The upper portion of the catchment corresponds to a high mountain landscape, with the highest concentration of small glaciers and ice fields of the mountain range [28].  Table 1 by their acronyms.
The estimated mean annual precipitation is 1840 mm [25], exceeding 2500 mm per year above 3000 m a.s.l. [29] Spring and autumn are the wettest seasons, while summer exhibits a marked dry period. Mean lapse rates have been estimated between -0.0049 and -0.0056 • C m −1 [30,31], and the 0 • C isotherm at around 1700 m a.s.l. during the cold season (November-April). This leads to a deep and extensive snow cover from late autumn to late spring. Consequently, the Ésera river exhibits a marked nival river regime with high flows in late spring and early summer. The lowest runoff occur in winter as consequence of snow accumulation, followed by late summer (August and September) in response to the lowest precipitation [32]. The annual runoff of the catchment is 309.2 hm 3 (9.8 m 3 s −1 on average).

Data
There are two pluviometers from the Automatic Hydrologic Information System of the Ebro river basin (SAIH) in the whole catchment. Temperature measurements are recorded at three snow gauges ("telenivometers"), which also record snow height and snow water equivalent (SWE). Two of the snow gauges are inside the catchment and a third one very close to the catchment boundary. The location of measuring stations is indicated in Table 1. Rain gauges are tipping bucket pluviometers, which are not suitable for solid precipitation and need careful calibration (WMO 2008). Snow height is recorded by an ultrasonic sensor and snow water equivalent using a cosmic-ray snow gauge. All data are public but additional information on sensor type and specifications is limited. These measured data were used to evaluate the model output, but the snow model was fed with WRF data directly. Precipitation distribution is difficult to evaluate as there were only two measurement sites, insufficient for a correct spatial extrapolation. Table 1. Coordinates of the measuring stations indicated in Figure 1. SWE is snow water equivalent, pcp is precipitation, temp is temperature and P is atmospheric pressure.

Station
Code Initial snow extent and albedo were derived from a Landsat image acquired from the U.S. Geological Survey and readily available through their online site http://earthexplorer.usgs.gov/. For the derivation of catchments and to calculate terrain parameters in the snow model, a Digital Elevation Model (DEM) was used. The original DEM at 5 m resolution was Lidar-generated and provided by the Spanish Geographical Office (IGN, Instituto Geográfico Nacional), this DEM was resampled at 30 m resolution in order to minimise artefacts, to make it compatible with the resolution of the Landsat images and to increase the computation speed of the snow model.
We approached the publicly funded Spanish Met Office (AEMET) for additional data, but were denied. We also tried to get data from the hydropower company operating the river dams, although this was an unsuccessful effort.

Remote Sensing
Initial albedo and snow cover extent were derived from a Landsat 8 image of 10 June, 2013. The snow extent was obtained by estimating the Normalised Difference Snow Index (NDSI), according to Equation 1 : where ρ n is TOA (Top Of the Atmosphere) reflectance of Landsat band number n [see for example : 33] For the albedo we used the Landsat surface reflectance and the conversion from narrow band to broadband was done according to Wang et al. [34]. Additional corrections were made for the angle of incidence of the sun on the slope surface, which can vary greatly on mountain areas. The solar-terrain geometry was derived from the DEM following Corripio [35] and using the R package insol [36]. The resulting initial albedo is shown in Figure 2 and compared the original false colour Landsat image. The northern fringes of the image are partially covered in clouds, in these areas the snow cover limit was set according to the elevation. The limit was derived from the surrounding cloud-free areas as the mean lower elevation with snow. This limit was dependent on the orientation of the slope, and was calculated in steps of 15 degrees from 0 to 345 degrees.

Snow Model
We use a multi-layered, distributed snow energy balance developed by Corripio [37], and applied successfully to previous studies in the Andes [38] and the Alps [39]. The energy balance of the snow pack can be described as in Equation 2: EB is the net energy balance at the snow pack, I ↓ is incoming short-wave radiation, α is snow albedo, L W↓ is received long-wave radiation, both from the atmosphere and surrounding slopes, L W↑ is emitted long-wave radiation, Q S is turbulent sensible heat transfer with the atmosphere, Q L is latent heat transfer, Q r is heat transfer due to precipitation, and Q g is subsurface heat transfer within the snow pack. If the snow pack reaches the melting point temperature EB will be the energy available for melting the snow, otherwise the net energy is used to change the snow pack temperature. Sublimation depends on the specific humidity gradient, whether the snow pack reaches melting point or not. Short wave incoming radiation was modified locally to account for shading, sky view factor and reflected radiation. A detailed description of the radiation components and the effect of surrounding topography is explained in Corripio [37]. Long Wave was derived from the WRF model, an additional cloud factor following Iziomon et al. [40] was applied to discriminate between the direct and diffuse radiation fraction depending on cloud cover. The WRF model calculates incoming long wave radiation independent of the local topography, so a further modification was applied to account for terrain influences such as sky view factor. Turbulent fluxes with the atmosphere were calculated according to the bulk aerodynamic flux method [41][42][43][44]. This method seems to be one of the most reliable according to a comparative study by Sexstone et al. [45]. Heat transfer due to precipitation was estimated according to Brun et al. [46]. The subsurface component of the model is a simplified three layer scheme where heat transfer is computed following Oke [47] with a simplified solar radiation snow extinction following Fukami et al. [48] and Warren [49]. Initial Albedo is derived from a Landsat image and local measurements, the temporal evolution of albedo is estimated according to an ageing parameterization considering local air temperature [50,51]. The temperature field was interpolated to the high resolution model according to local elevation and lapse rates from the WRF model. Snow-rain limit was prescribed by a wet bulb temperature of 1.3 • C [52,53]. This gives a 6 of 22 satisfactory result but an estimation that includes the precipitation rate could be an improvement. Finding the correct solution is not easy as the identification of the correct snow line is difficult during heavy precipitation and on mountain ranges with sparse instrumentation. The detailed equations of the model and derived variables are given in the Appendix.

Atmospheric Model
We run a nested WRF model ARW core [24], with three domains at 18, 6 and 2 km over the Pyrenees. The model initial and boundary conditions were derived from the GFS global model at 00Z every day and run for 30h to overlap the next run. The available resolution of the GFS at the time was 0.5 degrees, while the current resolution is 0.25 degrees. We did not consider using data from the European Center for Medium-Range Weather Forecast model (ECMWF) for the initial conditions as they are prohibitively expensive and beyond the budget available for this study. The nested model was run daily from the 11th to the 26th of June.
WRF model allows for different physical and dynamical settings and different radiative schemes, and these have an impact in the model ability to predict intense precipitation [54,55]. It would be desirable to make extensive tests to find which combination of settings perform better during heavy precipitation events in the Pyrenees. To our knowledge this is still to be tested. In the absence of this information, we decided to use the settings shown on table 2, based on information from previous work in the region and other mountain areas using the WRF model [56][57][58]. For a detailed description of the settings see the WRF User's Guide [59]. The output of the WRF model is in NetCDF format and the projection is in geographical longitude and latitude. To make it compatible with the DEM in Universal Transverse Mercator UTM projection, the files were reprojected using a Linux implementation of the GDAL -Geospatial Data Abstraction Library [60].

Discharge
From snowmelt and rainfall inputs river discharge was estimated using a parallel Single Linear Reservoir (SLR) model for every catchment. The SLR model assumes that storage S is linearly related to outflow Q by where K is the storage coefficient and has units of time [61]. Combining this equation with the hydrologic continuity equation [e.g.: 61,62] and integrating gives: where u is inflow at time t.
Hannah and Gurnell [63] suggest deriving the parameter K from the gradient of a semilogarithmic plot of discharge over time when there is no recharge, following Equation 5: But deriving the parameter K is problematic without accurate discharge measurements. The problem in this case is that the only river gauge is located after a dam used for hydroelectric generation. Thus, the river flow may change considerably depending on whether turbines are operating or not. If electricity is generated, the stored water in the dam is evacuated directly through a pipeline and bypass the river gauge. To complicate things further, there are two additional dams upriver that channel water to a power station located at the lower dam. The outflow of the lower damn for electricity generation is 30 m 3 s −1 ; the capacity of the lower damn (Linsoles) is 3.5 hm 3 ; the two upper dams (Estós and Paso Nuevo) are linked to the lower dam by a pipe with a maximum flow of of 36 m 3 s −1 . This situation makes difficult a proper calibration of the storage coefficient and therefore we have made an approximate calculation by trial and error. We used six different linear reservoirs corresponding to every sub-catchment, with different storage coefficients for snowmelt and rainfall. Travel times for every sub-catchment are taken into account and calculated with the GRASS module r.traveltime [64]. Singh et al. [2] recalls previous work showing that heavy rain may increase the speed of water drainage through snow, as more efficient channels develop in the snow pack. To account for this change we have modified the storage coefficient linearly with time.
The Aiguallut sub-basin on the north east of the main catchment (AG in Figure 1) is an endorheic basin, which drains to the Garona river, north of the area of study and already in French territory. It is a well developed karstic network with several intakes and an outflow at Uelhs deth Joeu to the Garona river [65,66]. Maximum measured outflow is 11-12 m 3 s −1 [65], at the time of the flood event the main intake at Aigualluts was saturated and discharge overflew to the Ésera river, thus this maximum value was subtracted form the main discharge of the sub-basin. Figure 3 shows a good agreement between measured temperature and temperature simulated by WRF, with an overall mean error of 0.21 • C and a mean absolute error of 1.51 • C at La Renclusa (point TN1 in reference table and map). However, maximum measured temperature is higher (up to 3.5 • C) than simulated temperature during days with high incoming solar radiation. The temperature probes are shielded but not ventilated, which may result in an overestimation of the maximum measured value [67]. The WRF model simulates the amount of precipitation approximately, as seen in Figure 4. It shows the maximum simulated precipitation to the west of maximum measured precipitation, but the sensor network is too sparse to decide whether there is a displacement bias or it is a realistic simulation. The simulated north -south gradient in precipitation is also observed in the measured data, between Llanos del Hospital (PV2) and Ampriu (PV1) The geomorphological evidence also suggest a strong gradient as the Vallibierna sub-catchment (VA) was showing less sediment transport and deposition than in previous floods events in the area (Santiago Somolinos, personal communication).   due snow accumulation is missed by the model. The Nash-Stuclife efficiency model coefficient [69] gives values ranging between 0.96 and 0.99, confirming the good agreement between measurements and simulations (Table reftab:03). Figure 5 also shows that except in the highest TN (Salenques at 2600 m a.s.l.) snow cover completely disappeared by the end of the period. Figure 6 shows the snow extent retrieved from a Landsat image in 26th June, 2013; and the simulated snow cover at the same date. It confirms the ability of the model to simulate the spatial distribution of snow cover depletion during the previous days and during the ROS event. Table 3. Nash-Sutcliffe efficiency values for the measured and simulated snow ablation at three snow gauges. During the day when peak flow occurred, the snow energy balance was dominated by a positive influx of energy from turbulent interchange with the atmosphere. Both sensible heat and latent heat transfer were positive, indicating condensation of atmospheric moisture on the snow (Figure 7). Shortwave radiation was second in importance, followed by advective heat transfer from the rain. Net long wave was positive, contrary to most days when there is an effective radiative cooling of the snow pack. Simulated subsurface fluxes were negligible during the ROS event, except for a short period early on the 19th of June when there was a positive flux toward the surface snow pack. For the calculation of heat transfer due to rain (Equation 33) it is necessary to know the temperature of the rain. Byers et al. [70] measured rain temperatures and found that after an intense precipitation event rain temperature was always close to air temperature (about 1 • C difference). For this study we assumed that the rain temperature was 1 • C colder than the air temperature. During the time of peak precipitation, simulated lapse rates above the surface were in the range between -0.003 and -0.0045 • C m −1 , thus one degree colder means that the rain had the temperature of a layer between 200 and 300 m above the surface.  Figure 8 shows the observed and simulated discharge at Eriste gauging station, the grey uncertainty shade is ±30m 3 s −1 corresponding to the intake of the hydroelectric plants. The peak discharge is well represented in the simulation (Nash-Sutcliffe coefficient of 0.76), but the river flow seems to decrease very slowly, with a secondary peak in precipitation that is larger than the observed variation. Rainfall is the main contribution to peak flow. Snow melt contribution is very well simulated before the ROS event, which adds confidence to the snow model results. Snowmelt was the main contributor to total water discharge during the previous days. Discharge was larger than the average for the season, probably causing saturation of soils and a very fast hydrological response to the rainfall event. During the flood event of the 18-19th of June estimated snowmelt represented 33% of total runoff, on the day of maximum precipitation (the 18th) snowmelt was calculated as 28% of total discharge and 23% at the time of peak flow.

Discussion
We believe that the combination of high resolution numerical weather prediction models, together with energy balance snow models are excellent tools to predict snowmelt and mountain river discharge, both during extreme events such as intense ROS, or during average conditions. We also think that is useful to stress uncertainty and limitations, so that we can improve this kind of simulation in the future. Simulated precipitation from the WRF model approximates the timing and intensity of the ROS event of 18-19th June, 2013 in the Ésera river. There are large discrepancies with local measurements ranging from 8 to 38% difference in accumulated precipitation. The sensor network in the area is too sparse to detect where the spatial distribution is correct or there is a displacement bias, as has been noted in other WRF simulations in the Pyrenees [58], as well as in other regions and for other variables [71]. Pennelly et al. [57] has shown the impact of cumulus parameterizations on the location and magnitude of precipitation events at lower resolution. It is unclear if this effect is propagated to the higher resolution nested domains, even when they can allow for resolving rather than parameterizing convection. The simulation of the precipitation peak is good enough for regional scale simulations, but it would be a great advantage if location accuracy could be improved. A thorough investigation into the best parameterization and best resolution for accurate numerical precipitation forecast in these mountain ranges is strongly recommended.
There is a strong north to south precipitation gradient, both in the measurements and the simulations. There is probably a vertical gradient in precipitation, too. This has been observed in many other mountain regions. For example Liston and Elder [72] use a seasonally variable linear gradient, where precipitation increases with altitude, according to work in the western US by Thornton et al. [73]. Closer to the Pyrenees, in the Alps Lehning et al. [74] have noted the problems of deriving a reliable precipitation gradient, either because of lack of measurements at high altitude [75,76], or because the difficulty of measuring solid precipitation at high elevation and with strong winds [77]. In fact, extensive work in the Swiss Alps reveals precipitation gradients ranging from -69 to +15 mm/km in summer or from -27 to +25 mm/km in winter [78]. We do not yet have enough information to derive a reliable elevation gradient of precipitation in the Pyrenees, and therefore it was not applied to the spatial interpolation of precipitation over the model domain.
The simulation of snow melt at a single point is successful, with a mean Nash-Sutcliffe efficiency value of 0.98 (Table 3) and suggests that the energy balance approach to simulate snow melt is a reliable one. This approach works well both for clear days, when solar radiation is the main energy input, and for overcast and rainy days, when the main energy input are turbulent fluxes with the atmosphere and there is an additional influx of long-wave thermal radiation. There is some discrepancy at the highest snow gauge where the model underestimates the amount of fresh snow on the 13th of June. This is probably due to a simplistic estimation of snow-rain limit which does not take into account the precipitation rate. WRF modelled temperature is very similar to measured temperature, but unfortunately we do not have measured data on short-wave radiation, long-wave radiation nor albedo. A full meteorological station measuring these variables would be of great help to test further the reliability of the WRF model. It is interesting to note that the mean slope of ablation rate changes very little after the 17th of June on the highest snow gauge, even when there is a drop of temperature of almost 10 degrees. The reason may be the increase in downwelling long-wave radiation from a cloudier atmosphere, which compensates for the reduced energy input from lower air temperature and reduced solar radiation. This highlights the remarks of Pomeroy et al. [79] on the unsuitability of simplistic temperature index models to simulate snow ablation and rain on snow episodes.
The distributed snowmelt model performs very well when comparing the total output of meltwater to discharge on days when there is no precipitation (Figure 7). This is reassuring and indicates that this modelling tools are also useful for standard operational river discharge forecasts in snow-covered regions. This is useful to optimise planning of hydropower production and to minimise conflicts in water usage. On days when there is heavy precipitation, the sum of rainfall and snow melt agrees well with measured discharge (Nash-Sutcliffe efficiency value of 0.76 overall), but rainfall decays too slowly on the WRF model. During heavy precipitation, rain creates more efficient channels and discharge is faster through snow [2]. An additional reason for this behaviour of the model could be an overestimation of the snow line altitude. Minder et al. [80] have reported observations where the snow line on the mountainside is lower than the estimated snow line in the free atmosphere and provide a theoretical explanation supported by model simulations. The main reasons for the depression of the snow line on the mountain slope are a) latent cooling from melting precipitation; b) microphysical melting distance (a 10 mm snow flake would melt 100 m lower than a 5 mm snowflake); and c) adiabatic cooling of orographically forced air masses. This situation has an important impact on runoff, for example in the Sierra Nevada of California a rise in the snow line of 610 m (1000 ft) would triple the runoff from three mountain catchments [81]. Ideally the storage coefficient to estimate discharge should change with time to reflect snow ripening and snow pack reduction, we have set it to vary linearly with time. Further calibration efforts are hindered by the artificial fluctuation of the river discharge in a catchment with hydropower use. This approach uses an empirical calibration, which is not ideal to simulate extreme events. For future work it would be advisable to get enough information on land cover and soil properties to implement different runoff models that are less sensitive to calibration (for example [79]).
At peak flow during the flood event, simulated discharge due to snow melt is 23% of total discharge, and 28% during the day of maximum precipitation, with rainfall being the main contribution to the flood. Problems caused in the catchment due to excess precipitation were exacerbated by a very high freezing level, which caused most of the precipitation to fall as rain, with only some snowfall towards the end of the event on the highest peaks. We did some tests with the GFS and the WRF models (not shown in this paper) to evaluate the lead times for a reliable warning in future cases. Both models identify well the height of the 0 • C isotherm up to 2.5 days before the event. GFS detects unusual precipitation but lacks resolution to identify the valleys most likely to be affected. WRF underestimates precipitation three days ahead and forecasts it approximately correct 24 h before the event.

Conclusions
We have shown that mesoscale numerical weather prediction models as the WRF are useful tools to predict the timing and location of rain on snow and flood events in mountain catchments of complex topography. Estimating the exact amount of precipitation is still the weakest point, although it can be improved by a systematic evaluation of current models spatial resolution and parameterization options [54]. An energy balance of the snow pack coupled with precise treatment of terrain parameters is able to simulate the snow depletion in the catchment and therefore the snowmelt input to discharge. Agreement between simulated and measured discharge is excellent in days without precipitation, which highlights the usefulness of these models for operational flow forecasting. The partition of energy fluxes show that the main input during this ROS event was turbulent fluxes of sensible and latent heat. Both fluxes were positive, implying strong condensation of atmospheric moisture. A simplified linear reservoir model to calculate the timing of the discharge gives satisfactory results, but can be improved by applying models that are less prone to calibration. This, however, requires extensive data on land cover, soil types and depth and snow depth. An additional source of uncertainty is the actual transition between snow and rain, we have used a simplified approach based on the local wet bulb temperature. This gives a satisfactory result but an estimation that includes the precipitation rate could be an improvement. Finding the correct solution is not easy as the identification of the correct snow line is difficult during heavy precipitation and on mountain ranges with sparse instrumentation. We have also found that both global and mesoscale models can give lead times of at least 2. which can be very helpful to to implement early warning systems and to avoid and mitigate potential flood damages.

Abbreviations
The

Appendix
Here we describe the detailed components of the energy balance of the snow pack. The incoming global solar radiation on every grid cell of the DEM is calculated as: where I n is direct normal radiation cos θ is the cosine of the angle of incidence of the sun on the terrain slope, f sh is a shading factor (0 if in shade, 1 if in sun), α is the snow albedo, I d is diffuse solar radiation, I dr is diffuse solar radiation reflected from the surrounding slopes, f skv is the sky view factor or fraction of upper hemisphere visible from the DEM cell, and f rsh is an estimation of the ratio of surrounding cells in the shadow. The ratio of direct to diffuse solar radiation is given by some compilations of the WRF or can be derived from a parametric solar radiation model, we used that of [82,83] and additional work by [84], [85], [86] for the estimation on cloudy days. The cosine of incidence of sun on the slope is the dot product s · n, where s is a unit vector in the direction of the sun: with ϕ as latitude, δ is declination, ω is the hour angle and n is the unit vector normal to the surface. For a regular DEM of cell size and elevation values z at row and column i, j respectively, n would be: For a detailed description of the sun terrain geometry and the calculation of sky view factor see Corripio [35]. Long wave downwelling radiation L ↓ is derived from the WRF model and compared to clear sky long-wave flux to determine the degree of cloud cover.
where is emissivity, T is temperature (K) with subscripts a and s for atmosphere and snow respectively and σ is the Stefan-Boltzman constant ( 5.6704x10 −8 Wm −2 K −4 ) From this, atmospheric emissivity is calculated according to Prata [87] as where precipitable water (w p ) is calculated following an empirical equation given also by Prata [87]: where T a is screen level air temperature in K. The actual vapour pressure is e 0 = e * RH, where e * is saturated vapour pressure computed following [88] polynomials (Equation 13), and RH is relative humidity from 0.0 to 1.0.
The profiles for wind speed, water vapour and heat, considering measurements at one level and at the surface, as: Subscripts m, v, h refer to momentum, water vapour and heat respectively. The Ψ s functions are stability corrections, which depend on the Obukhov's stability length. There is some disagreement on their values, which have been obtained empirically [89]. The values used here are similar to those used by Marks and Dozier [90], also reported by Brutsaert [89]. Thus, for unstable conditions (ζ = (z u − d 0 )/L ≤ 0): Ψ sh (ζ) = 2 ln 1 + χ 2 2 (20) with χ = (1 − 16ζ) 1/4 For stable conditions (ζ = (z − d 0 )/L > 0): where The transfer mechanisms of momentum and of other scalar admixtures are different at the surface, and consequently the roughness lengths have different values for momentum, water vapour and heat, which were calculated following Andreas [43]. He developed a fitting polynomial of the scalar roughness lengths over snow and sea ice (the coefficients for this are given in table 4): ln (z s /z 0 ) = b 0 + b 1 ln R e + b 2 (ln R e ) 2 (23) where z s refers to water vapour or heat depending on the coefficients used. The Reynolds number R e is a ratio of inertial to viscous forces [91] defined as: and ν = µ/ρ is the kinematic viscosity of air, with µ the coefficient of dynamic viscosity for air. Table 4. Values of the coefficients in the polynomials (equation: 23) that predict z 0 /z for temperature and water vapour [43] for smooth, transition and rough surfaces. The air density ρ for a given pressure is computed as the sum of densities of dry air ρ d and water vapour ρ v [89,92]: where the partial pressure of water vapour is p d = p z − e 0 and p z is atmospheric pressure, R d is the gas constant for dry air and 0.622 = 18.016/28.966 is the ratio of the molecular weights of water and dry air. The specific humidity at screen level q a and on the snow surface q s is calculated as: For a one dimensional system atmosphere -snow surface layer -subsurface snow pack the thermodynamic equation describing the conservation of energy, neglecting melting and refreezing, can be described as [93,94]: This equation is applied in a very simplified form following Oke [47] and implemented using finite differences: where ∆Q s is the flux density energy used to heat the snow at the surface, and C s is the heat capacity of the snow C s = ρ s c s , ρ s is snow density and c s is the specific heat of snow (2.09×10 3 J kg −1 K −1 ). As the main energy input comes from the surface, while the subsurface remains at a relatively stable temperature, we need to consider the thermal conductivity of the snow and the energy losses to subsurface layers: where κ Hs is the thermal diffusivity of the snow and k s its thermal conductivity. The volume of snow considered is that affected by the daily temperature changes at the surface. This temperature oscillation decreases exponentially with depth, and its amplitude at any given distance from the surface can be estimated according to Oke [47] by ∆T z = ∆ = T 0 e −z(πκ Hs P) 1/2 (32) where P is the wave period. Brutsaert [89] has shown that, roughly, 95% of the wave is damped at a depth of 3(2κ Hs /ω) 1/2 , where ω = 2π/P. Heat transfer due to precipitation Q r is estimated following Brun et al. [46] as: where ρ is rain water density, C p is specific heat of water, P r is precipitation rate and T s is temperature of the snow and T r temperature of the rain.