2.1.

H-SAF Dataset

As part of the EUMETSAT H-SAF, a near-real-time SWE dataset is produced and freely distributed.³⁴ The H-SAF SWE product has some similarities with the GlobSnow dataset³⁵ produced by the Finnish Meteorological Institute (FMI). However, a key feature of the H-SAF procedure is the separate treatment of flat/forest areas and mountain areas following the algorithms developed by the FMI and by the Turkish State Meteorological Services, respectively.³⁶ In summary, both institutes provide estimates of SWE for the entire European continent, but the optimized estimate of H-SAF is used for each of the two regions (flat/forest or mountain) according to a mountain map as determined using simple rules on mean altitude and standard deviation of the slope derived from the GTOPO30³⁷ digital elevation model.³⁸

The flat/forest areas component combines the sparse synoptic SD ground observations with passive microwave brightness temperature maps retrieved from the Special Sensor Microwave/Imager sensor on board the US Air Force Defense Meteorological Satellite Program series of satellites. The underling assimilation procedure³⁶ estimates the SWE by inverting a semiempirical snow emission model³⁹ and by combining these estimates with a weather station-based Kriging-interpolated SD map using a Bayesian spatial assimilation approach. A snow density value is applied to each grid cell to convert SD to SWE. Given that spaceborne microwave-based SWE estimates are reliable only for dry snow, the assimilation is not performed under wet snow conditions.⁴⁰

Over mountain areas, SWE estimates are based on spaceborne data only,³⁶ using an algorithm analogous to the one used for flat terrain. In addition, empirical estimates of grain size and density are considered more suitable over these terrains rather than ancillary information due to the common data scarcity. The mountain algorithm is designed for dry snow only; hence no SWE calculation is done if wet snow status is detected by the dedicated H-SAF product.

The merged product comprises daily maps in near-real-time (12 h time lag) covering Europe (25-75°N lat, 25°W-45°E lon) at a 0.25 deg spatial resolution, starting in 2013.³⁶ When ground observations are not available in time for a near-real-time delivery or when wet snow conditions are detected in mountain areas, the daily maps are masked over the corresponding pixels, resulting in partially missing data in the daily dataset. To maximize the data coverage at the monthly scale, monthly maps were obtained as maximum value composites for those months with at least three valid daily datasets available.

The dataset has been previously validated over Finland and Turkey.⁴⁰ Some key outcomes of the validation exercise of this product are an overall accuracy of monthly estimates [in terms of root-mean-square error (RMSE)] on the order of 30 to 40 mm and a clear tendency to underestimate SWE for values $>150\mathrm{mm}$ due to signal saturation. Also it is worth pointing out that most of the validations were performed against ground data collected between January and March to focus on dry snow only.

2.2.

CMC and ECMWF Datasets

The CMC provides a dataset of estimated monthly SWE⁴¹ starting in 1998 and regularly updated to the last full year (i.e., 2019 in this case).⁴² Since this dataset is produced only once per year, it can be considered a consolidated dataset over the northern hemisphere based on in situ observations that are not constrained by near-real-time requirements.

The monthly mean SWE values provided by CMC are used in this study. Those maps were obtained by applying to monthly mean SD maps the snow density values stored in a look-up table and derived from the Canadian snow-climate classification for different snow-climate classes.^43,44 The SD analysis is performed using a combination of in situ daily observations and optimal interpolation with a first-guess field generated from a simple snow accumulation and melt model based on temperature and forecasted precipitation.⁴⁵ In situ observations include data from synoptical stations and aviation reports.

Even if the look-up table values were derived from information acquired over Canada, the dataset includes monthly average maps for the entire canonical northern hemisphere snow season (namely October to June, the months when SD observations are available) on a standard 24-km polar stereographic grid. For this study, data were resampled over Europe on a lat/lon regular grid at a 0.25 deg spatial resolution using the nearest neighbor algorithm to obtain a dataset co-registered to the H-SAF grid.

A full quality assessment of the SD dataset at the base of this product is described in Ref. 45. Over Europe, the authors report an accuracy in the modeled SD of about 13 cm (RMSE), with a slightly negative bias of 5 cm.⁴⁵ Overall, the major errors of this product were found in Arctic regions (where no observations are available) and land areas with mean annual maximum snow accumulation exceeding 300 mm (e.g., Canadian arctic lands and Tibetan plateau).

In addition to the H-SAF and CMC datasets, the snow map provided by the ECMWF as part of the ERA5 reanalysis dataset was analyzed. Currently, the ERA5 dataset includes a large set of atmospheric, land, and oceanic climate variables covering the full globe with a grid of about 30 km on a reduced Gaussian grid (Tl639) and freely distributed within 5 days after real time. ERA5T data may still undergo a correction in a second step but no later than 3 months beyond real time.

In the ERA5 modeling framework, SWE is first simulated according to a physically based energy and mass-balance snow scheme.⁴⁶ The outputs of the snow analysis are then updated based on a two-dimensional optimal interpolation of station observations of SD and snow cover.^47,48

In this study, the monthly product distributed by ECMWF through the Copernicus Climate Data Store⁴⁹ is used. The raw data were cropped over Europe and regridded on the same $0.25\mathrm{lat}/\mathrm{lon}$ grid used by H-SAF. Monthly data between 1998 and 2019 were used in the analysis, constituted only by consolidated maps (i.e., no ERA5T data were used).

2.3.

Evaluation Strategy

Droughts are generally associated with anomalous conditions compared with the long-term climatology; hence all analyses in this study are focused on standardized $z$ scores (i.e., anomalies), defined as

The agreement of the $z$ score products derived from the H-SAF, CMC, and ERA5 data is evaluated by two analyses focused on: (1) snow detection over recurring snow-covered regions as defined by the CMC dataset and (2) consistency of the anomaly values of the three datasets. To assess the suitability of the H-SAF dataset for an implementation within the EDO operational system, threshold criteria are defined based on the indices derived during the two analyses.

The first analysis is made to detect the area that may be affected by snow drought events, i.e., areas with a systematic presence of at least a minimum amount of snow. To do that, a “snow month” is defined as when the monthly SWE is greater than the minimum threshold of 10 mm (analogous to the value used in H-SAF as an optimal accuracy⁴⁰). This minimum threshold allows for excluding regions with modeled SWE values that are $>0$ but still too small to be considered potentially affected by snow drought.

For each month, the number of past years with snow ($N_s$) are computed, and the grid cells where snow drought conditions can occur are obtained by defining a minimum number of years ($N_s>10$) in the full 1998 to 2019 CMC dataset. The CMC dataset is used here as baseline since it is the only one based on consolidated in situ data. The same procedure is then applied to both H-SAF and ERA5 datasets to perform a confusion matrix analysis, the outcomes of which are synthesized by the statistical indicators probability of detection (POD), probability of false detection (POFD), accuracy (ACC), and Heidke skill score (HSS), defined as

The POD and POFD indices have opposite behaviors, with the former being best at 1 and the latter at 0. For this study, we assume a POD of 0.8 as a threshold value to define a target consistency, representing an overlap on at least 80% of the areas detected by CMC. Analogously, a value of 0.2 is used as an acceptable threshold for POFD. In terms of accuracy and skill, target values of 0.8 and 0.5 are used as criteria for assessing the suitability of the product for inclusion in EDO.

In the second step, the consistency among the SWE anomalies modeled by H-SAF, CMC, and ERA5 is analyzed by means of a linear correlation analysis. In the absence of a clear SWE reference dataset for validation and given the uncertainty underling all three datasets, the correlation analysis is used to identify areas where significant inconsistencies among anomalies in the three datasets are detected.

Two main metrics are extracted from the correlation analysis: the Pearson correlation coefficient ($r$) and its statistical significance ($p$). It is worth mentioning that for standardized anomalies the $r$ coefficient represents both the correlation and the slope of the regression (the intercept is always zero); hence it also quantifies the accuracy and the bias. In particular, the percentage of the domain with positive correlation values and $p<0.05$ is used to identify the overall consistency between two models, and a value of 80% is assumed to be the minimum threshold for a satisfactory agreement.

In addition, the fraction of monthly values that fall outside a buffer region around the regression line is computed for each cell ($F_o$), assuming an amplitude of the buffer equal to $\pm 0.5$. The value of the buffer corresponds to the width of the main drought classes in commonly used standardized precipitation index classifications.⁵⁰ The criterion for suitability of the dataset is set to a maximum of $1/3$ of the data that can fall outside the two buffer lines.

3.1.

Analysis of the Spatial Coverage

The need for snow drought monitoring is limited to areas that have a temporal-consistent snow presence, i.e., where a deficit compared with the usual snow amount may constitute a concern for water managers. The temporal consistency of the CMC SWE grid data is summarized in Fig. 1, where the number of years ($N_s$) with $\mathrm{SWE}>10\mathrm{mm}$ is depicted for each month in the CMC snow season (between October and June). These maps show how the snow cover reaches its maximum around January/February and is negligible in October and June.

Fig. 1

Spatial distribution of the number of years ($N_s$) with $\mathrm{SWE}>10\mathrm{mm}$ for each month in the CMC snow season. The lower panel shows the temporal evolution of the percentage of the domain with at least 10 years with $\mathrm{SWE}>10\mathrm{mm}$ (out of a total of 22 years).

By defining a consistent snow presence when at least 50% of the years have $\mathrm{SWE}>10\mathrm{mm}$ ($N_s>10$, out of the total 22 years in the 1998 to 2019 CMC dataset), the temporal evolution of the percentage of the domain with consistent snow presence is reported in lower panel of Fig. 1, with values up to about 50% of Europe. The data show how two main areas are affected by consistent snow presence during the whole November to May period, namely the Alps and the Scandinavian peninsula, whereas the rest of the domain, mostly eastern Europe, have consistent snow only between December and April. Other threshold values (between 0 and 20 mm) were also tested (not reported) and confirm the general key outcome to focus the successive analyses on the period of November to May only (7 months).

To perform the same analysis on the H-SAF shorter time series, it is necessary to define a threshold analogous to 10 years (used for 1998 to 2019) for the period of 2013 to 2019 covered by this dataset. With this aim, threshold values between 2 and 5 years were tested on the CMC (not shown), selecting the value that returns the minimum difference in terms of year-average coverage. This optimal value, corresponding to $N_s>3\text{years}$, was used to estimate the overlap between CMC and the H-SAF and ERA5 datasets in the period of 2013 to 2019, as reported in Fig. 2(a).

Fig. 2

Overalp between CMC and both H-SAF (gray bars) and ERA5 (black bars) snow areas for the (a) whole domain and (b) mountain areas only.

These data show that both datasets detect about 85% of the grids defined by CMC as consistently snow-covered (7-month average of 84.9% and 85.3% for H-SAF and ERA5, respectively), with slightly more temporal-consistent results for ERA5 compared with H-SAF (7-month standard deviation of 5.8% and 4.0% for H-SAF and ERA5, respectively). The high overlap among all three datasets suggests an overall good agreement on the detection of the areas potentially affected by snow drought at the European scale, while the temporal consistency in the results further suggests a minimal impact of the selected threshold since lower SWE values may be expected early and late in the season due to SWE values being closer to the threshold.

One of the main features of the H-SAF algorithm is the coverage of mountain areas through a dedicated subalgorithm, the lack of which is a major limitation of other remotely sensed based products,²³ including the similar GlobSnow product. For this reason, the data in Fig. 2(b) report the overlaps in the case of mountain areas only, defined in this study as the common areas in the H-SAF and GlobSnow mountain masks³⁶ [see the mask map reported in Fig. 4(c)]. These data show a negligible deterioration of the overlap for the ERA5 dataset ($83.1\pm 4.8$) compared with the full domain, which is in contrast to the noticeable reduction observed for H-SAF ($36.6\pm 8.1$). These results seem to suggest that both CMC and ERA5 have quite consistent outcomes in terms of coverage, whereas the H-SAF product seems to be characterized by a noticeable incompleteness over snowy mountain regions, even if major snow areas are captured by the dataset in high-elevation regions. Part of the reason for the mismatch in the extension of the snow area in mountain regions may be ascribed to missing or low-quality satellite acquisitions (i.e., wet snow conditions, not modeled in mountain areas) rather than poorly detected snow conditions. This explanation is further supported by the worst performing months, November and May, which are located early and late in the season, respectively, when freshly accumulated wet snow is expected. Instead, it excludes a major effect of the selected 10 mm threshold since areas close to this limit are not predominantly located in mountain regions.

A more in-depth analysis of these SWE coverage maps was performed by means of a confusion matrix, the results of which are reported in Fig. 3 and summarized for the full study area in Table 1. Although the POD results are analogous to the overlap analysis reported in Fig. 2, the POFD highlights how a higher number of false detection is observed in the H-SAF product ($0.21\pm 0.05$) compared with ERA5 ($0.05\pm 0.02$), which seems to contrast with the previously reported underestimation of snow coverage over mountain regions, but it is in agreement with the small positive bias of SWE observed in H-SAF validation studies for low SWE values ($<150\mathrm{mm}$).⁴⁰ These results combine in a general better matching between CMC and ERA5 compared with H-SAF, with the latter showing fewer areas with consistent snow presence in mountain regions and a slightly higher extension of the area affected by consistent snow beyond the ones detected by CMC and ERA5.

Fig. 3

Indicators derived from the confusion matrix between CMC and both H-SAF (gray bars) and ERA5 (black bars) snow areas. Probability of detenction: (a) POD, the closer to 1 the better; (b) POFD, the closer to 0 the better; (c) ACC, the closer to 1 the better; and (d) HSS, positive values are skillful.

Table 1

Summary of the results of the confusion matrix analysis for the full study domain. See Fig. 3 for a synthetic description of the different indicators.

The main effect of these behaviors is nicely summarized by both ACC and HSS data, which are both higher for ERA5 ($0.92\pm 0.02$ and $0.78\pm 0.07$, respectively) compared with H-SAF ($0.82\pm 0.02$ and $0.55\pm 0.16$, respectively). Overall, all of the statistical metrics for the H-SAF product at the European scale are in line with the values set for a suitable use for operational drought monitoring, most notably the high values of POD and ACC (0.85 and 0.82, respectively), and they are generally comparable with the values obtained for ERA5. On average, the three models seem to be overall consistent in terms of the areas that can be considered consistently snow-covered, with the only exception being the still present underrepresentation of snow coverage in mountain areas in the H-SAF dataset.

3.2.

Analysis of the Consistency Among SWE Anomalies

Once the overall agreement between the snow-covered areas retrieved by CMC, H-SAF, and ERA5 products was established, the consistency of the three products in terms of SWE anomalies was evaluated by means of a correlation analysis on the period of 2013 to 2019 (up to 49 monthly values). Figure 4 reports the temporal Pearson correlation coefficient ($r$) maps for CMC versus H-SAF (a) and ERA5 (b), as well as ERA5 versus H-SAF (b). Each $r$ value was computed independently for every grid cell with at least 25 monthly values available.

Fig. 4

Spatial distribution of the Pearson correlation coefficients ($r$) between CMC and both (a) H-SAF and (b) ERA5 SWE anomalies, as well as between ERA5 and H-SAF (d). (c) The map depicts in blue the areas considered mountain terrain in the H-SAF algorithm.

The maps in Fig. 4 highlight an overall good correspondence among the SWE anomaly time series of the three datasets, which is particularly true over flat areas and eastern Europe. Sparse $r$ values lower than 0.2 can be observed over mountain regions [see image in Fig. 4(d)], including some negative values. On average, $r$ values are positive and statistically significant ($p<0.05$) in at least 90% of the domain in all cases (Table 2), highlighting an overall satisfactory consistency among the three SWE anomaly estimations over the majority of the domain. The lower average correlation observed between H-SAF and both CMC and ERA5 (compared with the correlation between CMC and ERA5) can be ascribed primarily to the reduced performances over mountain areas, as summarized by the average statistics over these regions (Table 3). The results over the mountain areas highlight a general low consistency among the outcomes of all three models, exemplified by both low average $r$ values and the large variability in the percentage of the domain with a positive and statistically significant correlation (ranging between 15% and 65%). It is worth highlighting that the impact of these inconsistencies on the overall comparison is downplayed by the limited extension of mountain areas.

Table 2

Summary of the correlation analysis between the three models’ SWE z scores.

Table 3

Summary of the correlation analysis between the three models’ SWE z scores over mountain areas.

Overall, a diminished accuracy over mountain regions may be expected for all three approaches. Although the difficulties that satellite products (e.g., H-SAF) have in reproducing SWE dynamics in mountain regions were been observed and discussed in Sec. 3.1, the accuracy of estimates from global reanalysis products (e.g., ERA5) also proved to be limited, with some studies demonstrating improvement only when regional nested approaches are implemented over specific regions.⁵¹ The reliability of the products based on the interpolation of in situ measurements (e.g., CMC), finally, is strongly linked to the availability of observations, which is generally completely lacking in high-mountain regions or only available once a year at the end of the accumulation period.⁵² Additionally, it has been shown that, in the case of high-resolution products, even interpolated data from snow pillows may not perform as well as energy-balance-based estimates over some regions.⁵³

The coarse spatial resolution of all three products may be another limiting factor over mountain terrains as it may not be suitable to capturing the complex spatial dynamic observed in such areas. Differences between point scale observation and grid-cell average SWE values have been shown to be up to 200% in regions with complex topography,⁵⁴ suggesting that the locations of in situ stations may not always be optimal for representing local average conditions.

Since CMC z scores are computed on a longer baseline period (1998 to 2019) compared with H-SAF and ERA5 (2013 to 2019), another potential source of discrepancy can be associated with this difference, so its impact was evaluated by computing the SWE anomalies for ERA5 over the period of 1998 to 2019 and by performing an additional correlation analysis on the full period (154 monthly values). The correlation map obtained in this case (not shown) is very similar to the one reported in Fig. 4(b), and the negligible differences in the results between the two periods are summarized by the synthetic data reported at the bottom of Table 2. Generally, the ERA5 average $r$ value is only marginally smaller when the full dataset is used, likely due to the increase in the sample size, but with a slightly larger percentage of the domain with the correlation being positive and statistically significant ($p<0.05$). This result seems to suggest a limited impact of the mismatch in the baseline periods.

An additional set of correlation analyses was performed to test the spatial consistency of the SWE $z$ scores, by computing the $r$ value for all grid cells in each of the 49 months in the period of 2013 to 2019. The average $r$ value for a specific month (from the $r$ values in the 7 different years) is reported in Fig. 5 for CMC versus H-SAF and ERA5 and for ERA5 versus H-SAF. These results confirm many similarities in the consistency between the three anomaly datasets. Even if the correlation between H-SAF and both ERA5 and CMC is lower than the $r$ value for CMC versus ERA5 (about 0.4 on average, compared with 0.6), the overall temporal dynamic is similar, with slightly lower values early/late in the snow season (i.e., December or May) compared with the central months. These low correlation values can be explained by the higher fraction of snow-covered cells in mountain areas during the early/late season compared with the middle season, as well as also by the difficulty of modeling wet snow for microwave-based approaches, which should be predominant during these periods. Expected lower performances for H-SAF during these months due to wet snow are also suggested by the product validation procedure, which focused on data collected between January and March only.⁴⁰

Fig. 5

Average spatial correlation for each month between CMC and both H-SAF (dark gray bars) and ERA5 (black bars) SWE anomalies, as well as between ERA5 and H-SAF (light gray).

The ending of the snow season, with a prevalence of snow melting, seems to be a particularly problematic period since minimum correlation values are observed in May in all three correlation analyses. During this stage, in addition to the already mentioned issue of wet snow modeling by satellite data, reanalysis approaches also are often limited by the large number of factors involved in the modeling of melting. Even if an improved scheme is implemented in ERA5 for a better timing of melting,⁴⁶ the improvement observed in the offline runs do not necessarily translate to the online outcomes.

A final test on the consistency of the three datasets was performed on the frequency, in which the SWE anomalies of H-SAF and ERA5 differ from the ones from CMC for more than the defined threshold of 0.5. Specifically, this fraction $F_o$ was computed independently for each cell [see example on Fig. 6(a)], and the frequency distribution across the domain is reported in Fig. 6(b) for both H-SAF and ERA5.

Fig. 6

Frequency analysis of data that fall outside the buffer region $\pm 0.5$ ($F_o$). (a) An example for a cell in Sweden, with the black line representing the linear regression and the dotted lines delimiting the buffer zone. (b) The frequency distribution for all cells in the domain, with the vertical dotted line highlighting the suitability threshold of 33%.

The analysis shows how high values of $F_o$ are relatively common, with about half of the domain with $F_o$ values greater than the threshold of 0.33 for H-SAF (i.e., only for half of the domain are the differences smaller than 0.5 for $<1/3$ of the monthly values). The agreement between CMC and ERA5 is slightly better (70% for $F_o=0.33$), and acceptable $F_o$ values can be observed over some mountain areas. By overlapping the areas with $F_o$ values below 0.33 for all three pairs of comparison, only about 30% of the domain has consistent outcomes in all three datasets according to the set criterion, mostly comprised of flat areas in eastern Europe. This result reinforces the difficulties in systematically achieving SWE anomalies that are consistent among all three products for the entire domain and further stresses the impossibility of achieving optimal results everywhere with a single dataset.⁵⁵