A virtual network for estimating daily new snow water equivalent and snow depth in the Swiss Alps

Abstract

Daily new snow water equivalent (HNW) and snow depth (HS) are of significant practical importance in cryospheric sciences such as snow hydrology and avalanche formation. In this study we present a virtual network (VN) for estimating HNW and HS on a regular mesh over Switzerland with a grid size of 7 km. The method is based on the HNW output data of the numerical weather prediction model COSMO-7, driving an external accumulation/melting routine. The verification of the VN shows that, on average, HNW can be estimated with a mean systematic bias close to 0 and an averaged absolute accuracy of 4.01 mm. The results are equivalent to the performance observed when comparing different automatic HNW point estimations with manual reference measurements. However, at the local scale, HS derived by the VN may significantly deviate from corresponding point measurements. We argue that the VN presented here may introduce promising cost-effective options as input for spatially distributed snow hydrological and avalanche risk management applications in the Swiss Alps.

1. Introduction

In cryospheric sciences such as snow hydrology, avalanche formation/dynamics and snow climatology, daily new snow water equivalent (HNW, mm) and snow depth (HS, m) are important measurement categories. Together with wind speed, HNW is of significant practical importance for avalanche hazard estimation (McClung and Schaerer, 1993; Egli, 2008) and represents a key input parameter for spatially distributed snow models (e.g. Lehning and others, 2006; Liston and Elder, 2006). The monitoring and modelling of water resources, represented by HS (respectively, snow water equivalent (SWE)), is the basic aim of snow hydrology (e.g. Grayson and Blöschl, 2001; Anderton and others, 2002). In particular, estimates of the spatial and temporal distribution of HS (or SWE) are essential (Luce and others, 1998; Skaugen, 2007) for water resources management (e.g. Schaefli and others, 2007) and flood prevention.

During the last decades, various expensive networks of manual and automatic point HS (and occasionally HNW) measurements have been set up in the Swiss Alps. Automatic methods have the advantage of providing data at high temporal resolution in terrain that is difficult to access. Durand and others (1993) presented a first analysis method for relevant meteorological parameters for snow models, while Egli and others (2009) investigated different methods to reveal a general feasibility benchmark in automatic estimation of HNW.

In this study, we present a virtual network (VN) to estimate HNW and HS on a regular grid (grid size 7 km) covering the entire Swiss Alps. The method is based on the HNW output data of the numerical weather prediction model (COSMO-7), developed by the Consortium for Small-Scale Modelling currently composed of the national weather services of Germany, Switzerland, Italy, Greece, Poland and Romania. The raw snow precipitation output of COSMO-7 drives a simple external accumulation/melting routine that allows HNW and HS to be computed for each gridcell. This routine does not account for individual physical processes (e.g. settling of snow cover and fresh snow). Instead, it evaluates bulk accumulation and melting rates based on COSMO-7 HNW data calibrated using observed HS data from existing snow-monitoring networks.

We believe that the VN presented here may introduce promising cost-effective options as input for spatially distributed snow hydrological models and avalanche risk management applications. If a numerical weather model and HS point measurement stations are available, a VN can be developed for other mountainous regions that are not as densely instrumented as the Swiss Alps.

2. Data

2.1. COSMO-7 model output data

A numerical weather prediction system is operated by MeteoSwiss (www.meteoswiss.ch) for a wide range of applications. This system is based on the COSMO-7 model (www.cosmo-model.org), which is a primitive equation model with non-hydrostatic, fully compressible dynamics. The prognostic variables include precipitation (separated into snow and rain) and many other meteorological parameters. COSMO-7 is used in two modes: (1) a free forecast mode, where the temporal evolution of the atmosphere is computed without any constraints other than the lateral boundary conditions from the driving model (European Centre for Medium-range Weather Forecasts, www.ecmwf.int/research), and (2) an assimilation mode to produce the best gridded representation of the atmosphere by blending available current observations with the model dynamic (the so-called data assimilation process). Note that no observed precipitation is assimilated in COSMO-7. The data assimilation process is not a simple interpolation of available observations, but instead produces a consistent state of the full atmosphere.

In this study we have chosen to use data from the assimilation mode of COSMO-7 between September 2005 and February 2008. For each day, snow precipitation (i.e. solid precipitation only) cumulated from 0000 UTC to 2400 UTC for the automatic reference stations (and from 0800 UTC to 0800 UTC on the following day for the manual reference stations) has been extracted from COSMO-7 and used to derive HNW_VN._RAW at the surface. Introduced abbreviations are summarized in Table 1.

Table 1. Table of nomenclature

2.2. Point measurements

For this study, 141 stations spread over the entire range of the Swiss Alps (Fig. 1) were available to provide daily measurements of HS during the period September 2005 to February 2008. The sites have been carefully selected in flat open terrain, with as little wind influence as possible, to ensure regional representativeness of the HS/HNW estimation (Egli, 2008). In order to derive HNW from HS measurements, we used a simple parameterization proposed by Egli and others (2009):

(1)

Fig. 1. Measurement stations: 141 point measurements of snow depths in the Alpine region over Switzerland. The stations are located at altitudes 800–3130 ma.s.l., where the three different elevation zones are indicated with squares, points and triangles. Reference stations are marked with a cross.

where δHS_MEAS._24h denotes the 24 hour difference between an HS reading (measured here in centimetres) and the respective observations from the previous day. Possible limitations of using Equation (1) are discussed in section 4.1. Note that for automatic measurements from the network of the Intercantonal Measurement and Information System (IMIS; see Rhyner and others, 2002), HS estimations were provided at 0000 h while manual HS readings (for observations, see Marty, 2008) were provided at around 0800 h. The values of HS were later checked for plausibility, and single missing/faulty values were corrected manually by interpolation. Finally, HNW was only calculated for the period 1 November–30 April for each season (analogous to Egli and others, 2009), while HS is considered during the entire seasons of 2005–08 if HS >0.

3. Methodology

3.1. Virtual network (VN)

HNW estimated by VN

The cumulated snow precipitation output of COSMO-7 (HNW_VN._RAW) constituted the basis of the VN. Thirty-five snow stations (of 141) were set aside for adjusting HNW_VN._RAW to measured HS (HNW_MEAS). These stations are referred to as reference stations (Fig. 1, crosses), while the remaining 106 stations are referred to as control stations. All stations were separated into three different elevation bands (800–1500 m, 1501–2200 m and 2201–3200 m; see legend of Fig. 1).

Reference stations were selected to cover all regions and elevation bands approximately evenly, in the horizontal as well as in the vertical space. COSMO-7 mesh points were assigned to corresponding snow stations by searching the gridpoint with centre nearest to the station. The horizontal station coordinates deviated about ±2.7km on average (maximum from the centre of HNWVN.RAW, while the vertical difference between stations and respective HNW_{V N.RAW} gridcells deviated by about ±250 m (Schaub, 2007). The VN therefore covers an area of approximately 7 km × 7 km on the surface (horizontal). In comparison with the point measurement stations, a vertical distension of about 5 0 0 m is provided.

The correction routine for HNW_VN.RAW first identified possible outliers of HNW_VN.RAW if HNW_VN.RAW exceeded the maximum of HNWr_MEAS at the reference stations (index r). Outliers were replaced by the maximum of HNW^r _MEAS for eachelevation band. Secondly, a correction matrix CorrMat_r was calculated from the difference between HNW_MEAS and HNW_VN._RAW for each reference station and for each elevation band:

(2)

For all remaining locations (index c) of the respective elevation band, the values of the correction matrix were then interpolated by an inverse distance interpolation:

(3)

where

(4)

is the inverse distance weighting function. d represents the Euclidean metric distance operator and N is the total number of reference locations (N = 35 here). The corrected COSMO-7 HNW estimation for the control stations (HNW^cVN _. CORR₎ was then defined by:

(5)

where negative values of HNWcVN.CORR were set to zero.

HS estimation by the VN

HS is basically calculated from the summation of HNW during the accumulation period, while a simple melting procedure is used to calculate HS during the ablation period. To calculate a temporary term TEMP.HS_cV_rN(t) for the COSMO-7 gridpoints, the values of HNW_cV_rN._CORR were added iteratively day by day:

(6)

Note that instead of a simple cumulation of HNW_VcrN._CORR(t) for each day, this iterative procedure also allows HS to be determined for snowfall events during the melting period. However, at this point TEMP.HS_cVrN(t) does not indicate a real snow depth measure. In order to convert TEMP.HS_cVrN(t) to corresponding HS data, we employed data from the reference stations to calibrate TEMP.HS_cVrN(t) using a simple linear parameterization:

(7)

Again, such a parameterization was obtained for each elevation band separately. The resulting parameters of Equation (7) (slope and intercept) were subsequently used to convert TEMP.HS to HS at the control stations:

(8)

The above procedure was applied for all days (considered here as the period of accumulation) except during ablation, which is defined as all days if mean(HS_rMEAS(t)) - mean(HS_rMEAS(t - 4)) < 0 and mean(HS^rMEAS(t)) < 0.75 × max(mean(HS^r _MEAS(t))). Note that this criterion does not take into account any information about the physical processes of melting (e.g. liquid water content). It approximately defines the point where the mean snow depth at the reference stations for each elevation band decreases due to melting conditions (Egli and Jonas, 2009). In case of melting conditions, we calculated a melting matrix from HS_rMEAS(t) - HS^r _MEAS(t - 1) for every elevation band at the reference stations:

(9)

We calculated melt rates for all remaining locations using the inverse distance weighting method, in the same way as for Equations (3) and (4):

(10)

where

(11)

Finally, HS (for melting conditions) was derived for the control stations as

(12)

Note that at the very end of the winter season the reference stations would eventually melt out, preventing any melting rates from being derived. In this case, the final melting rates before meltout have been used for the interpolation. In addition, negative values of HSVN have been set to 0.

3.2. Parameters of verification

HNW verification

In order to validate the performance of HNWVN at the control stations (in the following the index c is not included) by comparing to the point measurements HNWMEAS, and to compare the results to the performance of other automatic methods, parameters of comparison were used (Egli and Jonas, 2009). The parameters are briefly summarized below.

1. 1. Systematic bias :

(13)

1. 2. The absolute accuracy σ(δHNW):

(14)

1. 3. : The coefficient of correlation between log-transformed HNWV N and HNWMEAS.

2. 4. POD and FAR: Parameters adopted from severe weather forecast theory (Murphy and Winkler, 1987) have been used to detect certain classes of snowfall events, namely the probability of detection (POD) and the false alarm rate (FAR). To calculate POD and FAR, we used contingency tables (Table 2) for four classes of HNW intensity (Snow-NoSnow: HNWMEAS ≥ 0 mm and < 1 mm; Low: HNWMEAS ≥ 1mm and < 15 mm; Medium: HNWMEAS ≥ 15 mm and < 30 mm; High: HNWMEAS ≥ 30 mm).

   The POD and FAR were calculated for every class using

   (15)

   

   and

   (16)

   

   where b, c and d were derived from the contingency table (Table 2).

3. 5. Ranking points: A ranking point scale for the overall assessment relative to other HNW estimation methods was tested in Egli and others (2009). For each comparative measure as described above, the best performance was rated with 7 points and the lowest performance with 0 points. The different ranking points for each parameter and class were added.

Table 2. Contingency table for the validation of the HNW estimations of the virtual network. The point measurements at the control sites HNWMEAS are observed and the estimations of the virtual network HNWVN are forecast

HS verification

HS data from the VN (HS_VN) were validated against HS_MEAS by applying the following parameters of comparison (Egli and Jonas, 2009), restricted to the control stations.

1. 1. Systematic bias (α): The percent systematic bias α was determined by a least-squares fit to

   (17)

   

   where the standard error of α was also determined (errorα). Here, errorα is a measure of the spread of the data around the fitting line.

2. 2. Absolute error (RMSE): The averaged absolute difference between HSVN and HSMEAS, the root-mean-squared error (RMSE), was included in the comparison statistics:

   (18)

   

   where i indexes one day of the analysed periods if either HSVNi > 0 or HSMEASi > 0. Since HS is investigated instead of SWE (Egli and others, 2009), a direct comparison of these studies by a ranking point score is not appropriate.

4. Results and Discussion

4.1. HNW

Results of the validation of HNWVN and HNWMEAS are listed in Table 3 and presented in Figure 2. The parameters of comparison were calculated for both HNWVN.CORR and HNWVN.RAW compared to HNWMEAS. For reference, they are listed together with corresponding results of other HNW estimation methods which were tested in Egli and others (2009). The parameters for VN.CORR and VN.RAW are averaged over all 106 control stations (Table 3: nations), where the lines below and above the end of the vertical bars in Figure 2 (as well as the numbers in brackets in Table 3) represent ±half of the standard deviation.

Table 3. Comparative statistics to assess the performance of VN relative to competing methods tested in Egli and others (2009); notation of the systematic bias the absolute accuracy (σ(5HNW)), the coefficient of correlation and the ranking point system are described in section 3.2

Fig. 2. POD and FAR values for the four different classes of HNW intensity (diamond, point, triangle, star). The values are listed either for single point measurements or the average of 106 control stations, where the error bars denote the deviation from the mean.

Table 3 summarizes the parameters for HNW comparison ( σ(δHNW) and ) in order of the highest ranking points. While indicates a general underestimation of HNW _VN to HNWMEAS, denote a statistical overestimation. The results for VN.CORR and VN.RAW reveal a systematic bias close to 0. However, the considerably large deviation around the mean values of VN.CORR and VN.RAW indicates that some stations exhibit a large under-/overestimation of HNWVN.

A more detailed analysis of the frequency distribution showed that the individual at each control station are approximately equally distributed around the mean value of -0.03 mm. This indicates an arbitrary character which stems from the comparison of the point measurements HNWMEAS with data representing a much larger domain of estimation for HNWVN.CORR or HNWVN.RAW. This is supported by the fact that the stations are randomly distributed in a cuboid of 7 km × 7 km by 500 m (Schaub, 2007), which also explains why the deviations are generally similar for VN.CORR and VN.RAW.

For VN.CORR, σ(δHNW) demonstrates considerably small values for the mean (σ(δHNW) = 4.01 mm) with a deviation from the mean of about ±0.62. This low absolute error σ(δHNW) implies that, in principle, HNW_VN CORR can also be calibrated for each single station of the measuring network. The systematic bias of each station can therefore be reduced, resulting in a lower deviation (±) in . However, due to the comparison of a point measurement with a cuboid area, a calibration of HNWVN.CORR is not meaningful. In fact, the values of HNWVN.CORR may be considered as an average measurement of HNW over a large area, therefore representing a complementary network to the point measurements.

Taking the performance of SNOWPACK for σ(δHNW) as a benchmark, about 82% of the 106 control stations of the VN exhibit a σ(δHNW) < 5.07mm. The absolute accuracy between HNWMEAS and HNWVN.CORR is therefore comparable to the operational SNOWPACK-HNW estimations for manual reference measurements investigated at a single location. Note that SNOWPACK calculations are in operational use for Swiss avalanche warning (Rhyner and others, 2002) and provide different types of snowpack properties (Lehning and others, 2002). However, they are restricted to the locations of the IMIS network (Rhyner and others, 2002). In contrast, VN.CORR constitutes a complementary network over the entire Swiss Alps (although restricted solely to the HNW (and HS) estimation).

The POD/FAR statistics show that VN.CORR performs for all classes of intensities with better POD/FAR values than VN.RAW. Note that POD denotes the percentage of HNWMEAS, also measured by HNWVN, while FAR represents the percentage of HNW_VN which are not observed by HNW_MEAS. VN.CORR has similar features to the POD/FAR statistics of SNOWPACK, in particular for the two highest classes of intensity (Medium and High) and for the Snow-NoSnow class, which are the most important for avalanche risk management decisions. The good performance of VN.CORR in POD/FAR statistics with respect to the Snow-NoSnow class also implies that COSMO-7 is capable of distinguishing between snow and rain. Note that results do not deteriorate if evaluated for the three elevation bands separately (data not shown).

Furthermore, the calibration routine to derive VN.CORR implies that at the reference station HNW_VN._CORR = HNW_MEAS, independent of the processes leading to HNW_MEAS such as combinations of solid, liquid precipitation and evaporation.

The low performance of VN.CORR in POD (0.45) and FAR (0.35) for the Low class diverges from the SNOWPACK results. In this class, VN.CORR is comparable to the performance of the SIMPLE-HNW method. This is unsurprising given that HNW are determined using SIMPLE-HNW at the control stations. This limits the potential of our approach to deal with situations in which the change in snow depth is strongly influenced by both snowfall and settling. While Equation (1) implicitly includes the effect of settling of new snow and the entire snow cover, it may become inaccurate after the first day of a multi-day snowfall event. In such situations VN.CORR therefore may not perform as well as a physical snowpack model. Note also that the difference in the time of measurements between the manual and automatic measurements (8 hours; see section 2) may lead to erroneous HNW estimates in a few cases (e.g. if a significant precipitation event occurs between 0000 and 0800 h). However, an HNW analysis of manual and automatic stations separately (data not shown) revealed very similar results with respect to POD/FAR statistics, which implies that this effect is minor.

Overall, VN.CORR demonstrates an advantage over VN.RAW. In particular, the calculation described by Equations (1–5) demonstrated the largest improvement. The method of removing the outliers has a minor impact on the results of VN.CORR. Moreover, in comparison to the other methods, VN.CORR yields the highest number of ranking points and therefore constitutes an adequate automatic estimation for HNW over the entire Swiss Alps.

4.2. HS

Figure 3 shows the temporal development of the modelled HS (HS_VN, black curves) and the measured HS (HS_MEAS, grey curves) during winter 2006/07 for three exemplary locations of automatic measurement stations (ELM 2, DAV 2 and DAV 3). A temporal development of HS is first displayed where HS_VN describes a curve consistently below HS_MEAS (ELM 2), resulting in a systematic underestimation of α = 0.65 and a mean absolute error (RMSE) of 0.42 m. Secondly, a trajectory is shown (DAV 3) where the systematic bias (α) is close to 1 and RMSE = 0.19 m, indicating a strong congruence between the curves. Finally, a station (DAV 2) is shown which describes a path constantly above HS_MEAS, resulting in α = 1.54 and RMSE = 0.45 m (i.e. similar to that of ELM 2). Despite the fact that HS_VN shows considerable deviation from HS_MEAS for some stations, all stations generally reproduce the temporal progression qualitatively well. This is particularly true for the peaks during accumulation and the following settling (produced by Equations (7) and (8)). The same result is obtained for the ablation derived by Equations (9–12).

Fig. 3. The temporal development of snow depths derived by the virtual network (black curves) and the point measurements (grey curves) for three selected stations. The systematic bias (α) and the absolute error (RMSE) are indicated in the legend.

Basically, the ablation period is generated by the VN only by means of differences of HS, where the melt rate matrix (Equation (8)) actually refers to a combination of some fresh snow accumulation, settling and melting. As the peaks are attributed to snow accumulation by daily new snowfall (HNWVN.CORR), the deviation from HSMEAS stems from possible under-/overestimations of HNWVN.CORR, as discussed in section 4.1. Analogous to the results and discussion of the systematic bias of HNWVN.CORR, the specific values of a for each station are approximately equally distributed around the mean value of . Accordingly, while α is reasonably close to 1, and are very large (see Table 4).

Table 4. The parameters of HSVN evaluation statistics (α, errorα and RMSE) as described in section 3.2 calculated either for the analysis per station averaged over 106 stations (HSVN per station), or for the analysis of HSVN averaged over all stations (average of HSVN)

4.3. Spatial analysis

Since the analysis of the single point measurements of HSVN is limited due to its random nature, an investigation of HSVN and HNWVN.CORR regarding the Swiss Alps as an entire area is also considered. First, the mean of all snow depths of point and VN estimations ( and ) and their standard deviation (CT(HSMEAS) and σ(HSVN)) are calculated. Figure 4 displays the evolution of σ(HSVN) with (black symbols) and σ(HSMEAS) with (grey symbols), where each point of the trajectory represents one day of the year 2006/07. The curves show a characteristic hysteretic dynamics as discussed in Egli and Jonas (2009). The trajectories of the period of accumulation (points) and ablation (triangles) are clearly separated. The quasi-linear increase of σ(HSVN) with increasing highlights that accumulation of snow leads to an increase in the differences between sites. The path during ablation, on the other hand, is mainly attributed to the spatial distribution of melting rates in the Alpine region.

Fig. 4. The evolution of the mean snow depths (HS) and the standard deviation of snow depths derived by the virtual network (HSVN, black symbols) and the point measurements (HS_MEAS, grey symbols). The trajectory during the period of accumulation is indicated by points; the trajectory during ablation is indicated by triangles.

It has been shown (Egli and Jonas, 2009) that the curves of the hysteretic dynamics are similar between years and therefore characterize the seasonal development of HS in the Swiss Alps. Figure 4 shows that both the periods of accumulation and ablation are well reproduced by HSVN when compared to HSMEAS. The seasonal characteristic of the HS development is therefore also displayed by the VN. We therefore speculate that HSVN correctly reproduces the total amount of snow depths during a season at Swiss Alpine scale. Moreover, when evaluating a, errorα and RMSE for and (HS is averaged day by day over all stations), the percentage bias is as small as -4% (Table 4).

Furthermore, the variogram analysis of HNW (Egli, 2008) showed that the statistical correlation length of HNW over the entire Swiss Alps is about 50–60 km. The same correlation length was found in Schaub (2007), analysing the raw HNW output of COSMO (HNWVN.RAW). This additionally supports the assumption that HNWVN reproduces the regional precipitation patterns of the Swiss Alps appropriately and can be used as an input parameter for spatially distributed snow models.

Finally, the POD statistics of point HNW measurements as a function of the distance between two measuring sites (Egli, 2008) showed that, for the smallest scale (5–10 km), the POD is about 60% for HNW intensities more than 30 mm (High). The same value has been derived from a comparison of HNWVN.CORR and HNWMEAS. Again, a better performance of the VN HNW estimation in the POD statistics of two point measurement stations about 7 km distant would have been astonishing, since the large area of HNWVN.CORR covers a gridcell of 7 km × 7 km. As the POD decreases rapidly with the distance from a point measurement station (Egli, 2008), the VN may also be applied for avalanche risk management warning for the locations between the point measurements where no HNW estimation is available.

5. Conclusion and Outlook

In this study, the HNW output of the numerical weather prediction model COSMO-7 was coupled with a simple snow accumulation/melting model in order to provide HNW/HS grids of 7 km resolution for the entire Swiss Alps. The results for the HNWVN estimation showed that, on average, its performance is comparable to the performance of different automatic point measurements. HNWVN therefore represents a complementary network for avalanche risk management applications and can also be used as input data for spatially distributed snow hydrological models where HNW is required.

This is also the case for the HS estimation by the VN, where the statistical dynamics of the mean of HS and its standard deviation (σ(HSVN)) are congruent to the measured dynamics. As a consequence, the estimations of HS by the VN over the entire Swiss Alps may be used to estimate the total amount of snow and snow water equivalent for larger catchments. At a local scale, however, a direct comparison of the point measurement to the corresponding gridcell of the VN may result in considerable deviation. This may stem from the fact that a measurement over a large area (7 km × 7 km for the VN) is compared to a single point measurement.

Future effort is necessary to investigate the spatial scales over which the VN is capable of providing a good performance for estimation of HNW/HS. For this purpose we will extend the application to COSMO-2, a new high-resolution version of COSMO available from February 2008 for a mesh size of 2.2 km. Additionally, the differences in the VN between the assimilation mode and the forecast mode of COSMO-2 will be investigated. Finally, the principle of the VN presented here can be applied to other regions which are less densely equipped with automatic or manual point measurement stations than the Swiss Alps. If point measurement stations and a numerical weather prediction model are available, HNW and HS can be estimated with a refined accuracy for spatially larger extended regions for avalanche-risk and water-resources management.