The Crocus snowpack model within the Interactions between Soil–Biosphere–Atmosphere (ISBA) land surface model was run over northern Eurasia from 1979 to 1993, using forcing data extracted from hydrometeorological datasets and meteorological reanalyses. Simulated snow depth, snow water equivalent, and density over open fields were compared with local observations from over 1000 monitoring sites, available either once a day or three times per month. The best performance is obtained with European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Re-Analysis (ERA-Interim). Provided blowing snow sublimation is taken into account, the simulations show a small bias and high correlations in terms of snow depth, snow water equivalent, and density. Local snow cover durations as well as the onset and vanishing dates of continuous snow cover are also well reproduced. A major result is that the overall performance of the simulations is very similar to the performance of existing gridded snow products, which, in contrast, assimilate local snow depth observations. Soil temperature at 20-cm depth is reasonably well simulated. The methodology developed in this study is an efficient way to evaluate different meteorological datasets, especially in terms of snow precipitation. It reveals that the temporal disaggregation of monthly precipitation in the hydrometeorological dataset from Princeton University significantly impacts the rain–snow partitioning, deteriorating the simulation of the onset of snow cover as well as snow depth throughout the cold season.
Scientific interest in snow cover dynamics is continuously growing, particularly in relation to climate variability. Indeed, many physical processes controlling the climate system are strongly affected by the presence of snow on the ground. Because of its unique physical properties (high albedo, high emissivity, low thermal conductivity, latent heat sink while melting, and low roughness), the snowpack drastically changes the surface radiative and turbulent fluxes.
In the Northern Hemisphere, it strongly influences the land surface energy budget and thus the state of the lower troposphere (Cohen and Rind 1991; Groisman et al. 1994). Moreover, as Northern Hemisphere snow extent is characterized by high interannual variability (Frei and Robinson 1999), it has a potential predictive value for seasonal forecasting. Some studies have explored the potential of snow cover extent as a source of predictability in global circulation models (GCMs) (Schlosser and Mocko 2003; Peings et al. 2010). They highlight the importance of initialization in terms of snow cover extent for accurately reproducing the springtime temperature anomalies over Europe and North America. In addition to its local impact, snow is expected to influence remote large-scale atmospheric modes of variability. In particular, several studies suggest that Siberian snow cover extent in autumn is linked to the subsequent Arctic Oscillation phase in winter. This teleconnection is supported both by empirical studies (Cohen and Entekhabi 1999; Bojariu and Gimeno 2003) and by numerical sensitivity experiments (Gong et al. 2003; Fletcher et al. 2009; Peings et al. 2012). Other studies suggest a relationship between the Siberian snow cover extent and the North Pacific atmospheric circulation (Walsh and Ross 1988; Walland and Simmonds 1996; Orsolini and Kvamsto 2009). The relationship between the Eurasian snow cover and the Asian/Indian summer monsoon has been extensively studied as well (Hahn and Shukla 1976; Barnett et al. 1989; Bamzai and Shukla 1999). However, it is still a matter of debate, as the teleconnection is not stationary in the observed record (Robock et al. 2003; Peings and Douville 2009). An accurate representation of snow in the GCMs is therefore important in order to reproduce some large-scale atmospheric modes of variability in a consistent manner.
Beyond its impact on climate variability, the Eurasian snow cover is a major component of the global water cycle. Indeed, the direct runoff of snowmelt over frozen soil is an important contributor to the river streamflow in Siberia (Poutou et al. 2004; Grippa et al. 2005; Niu and Yang 2006; Decharme and Douville 2007). The snow cover partially insulates the soil surface from the atmosphere, thereby impacting the seasonal freezing depth (Yang et al. 2002; Zhang et al. 2005; Sushama et al. 2006).
The snow cover is also a major source of positive feedback in the climate system. However, most GCMs still suffer from recognized deficiencies in the simulation of snow, including its impact on Earth’s radiative budget (Roesch 2006). This in turn contributes to uncertainties in climate change scenarios. Therefore, many efforts have recently been dedicated to a better representation of the snowpack, mainly based on the implementation of multilayer physically based snow schemes in GCM land surface models (Oleson et al. 2010; Walters et al. 2011). The earth system model developed by the Centre National de Recherches Météorologiques (CNRM) and the Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique (CERFACS), referred to as CNRM-CM (Voldoire et al. 2013) and used for the Climate Model Intercomparison Project 5 (CMIP5), is still using a rather simple one-layer snow scheme, though more sophisticated multilayer snow models are currently being developed. As a first step toward the implementation of a more detailed snowpack scheme in CNRM-CM, the present study evaluates the performance of snow cover simulations over northern Eurasia driven by hydrometeorological datasets and meteorological reanalysis. It uses the Crocus snowpack scheme coupled with the Interactions between Soil–Biosphere–Atmosphere (ISBA) land surface scheme used in CNRM-CM (Noilhan and Mahfouf 1996; Voldoire et al. 2013; Vionnet et al. 2012). Recent studies have performed detailed snowcover simulations over northern Eurasia. Yang et al. (2010) have highlighted the impact of latitude and blowing snow parameterization on snow sublimation. Troy et al. (2011) have demonstrated the capacity of the Variable Infiltration Capacity (VIC) land surface model (LSM) to reproduce monthly averaged ground-based snow depth observations and analyzed the impact of different forcing datasets on the simulated water budget. Liston and Hiemstra (2011) used gridded weekly and monthly snow products over the last three decades to evaluate a high-resolution distributed snow model over a pan-Arctic domain. In the present study, a major difference lies in the detailed comparison between the simulations and the observed station datasets. Special attention has been paid to the capacity of the model to reproduce the observed daily snow depth series and the observed 10-day series of snow water equivalent and density over northern Eurasia. Several hundred individual stations from the former Soviet Union, extracted from the meteorological synoptic and hydrological survey observation networks, have been considered. The present study is limited to snowpack simulations over open fields, for which a considerable amount of ground-based observations are available. The sensitivity of the results to the driving dataset has been analyzed, leading to the identification of weaknesses in some meteorological variables.
2. Materials and methods
a. General methodology
This study focuses on the capacity of the Crocus physically based snowpack model, coupled with the ISBA LSM, to simulate several important characteristics of the snowpack [depth, snow water equivalent (SWE), density, duration, onset, and vanishing dates of continuous ground coverage]. The domain covers the watersheds of several major northern Eurasian rivers (Lena, Yenissei, Ob, Volga, and Amur) and the study period extends from 1979 to 1993. The snow simulations were evaluated against snow depth, SWE, and density observations from all stations within the considered domain available in the Historical Soviet Daily Snow Depth (HSDSD) and Former Soviet Union Hydrological Snow Surveys (FSUHSS) datasets. The meteorological variables necessary to drive the LSM were extracted from global hydrometeorological datasets and reanalysis. A correction was applied to the original meteorological data in order to account for the difference in elevation between the stations and their nearest grid points. Then different configurations of the LSM (different forcings and different Crocus options) were run over every station and compared with the corresponding observations. We focused on the bias, correlation, and root-mean-square errors (RMSE) of the snowpack characteristics mentioned above. In addition to the local simulations, we used one of the configurations to perform a 2D simulation at 0.5° resolution over the whole domain. The latter simulation was used to illustrate the spatial variability of the simulated snowpack.
b. Snow cover observations
There are many sources of observations related to snow cover over northern Eurasia. Several global products derived from satellite observations are available, but they are not very suitable for evaluating daily local snow simulations; for example, snow extent products based on sensors operating in the visible range of the solar spectrum are limited to indicating the presence of snow on the ground and are affected by cloudiness. SWE products, based primarily on microwave observations, show relatively large uncertainties, especially during the melting period (Pulliainen 2006) and in the case of deep snowpacks (Takala et al. 2011). They are also affected by the presence of vegetation in the considered pixels (Pulliainen 2006). Hence we focused on ground snow observations carried out at local stations.
The first source of snow observations comes from the HSDSD version 2 dataset (Armstrong 2001). It contains data collected in the former Soviet Union and processed by the State Hydrometeorological Service (Obninsk, Russia). They are accessible through the National Snow and Ice Data Center (NSIDC) data portal. The series include daily snow depth data measured at synoptic stations following the World Meteorological Organization (WMO) standards. Consequently, they were collected in open areas or clearings even in forest regions. WMO requires regular grass cutting, which means that snow depth should be measured over ground covered by very short grass. This is an important point that will be discussed further below. This dataset has two major advantages: its quality control and the fact that snow depth is recorded year round, even during the warm season when snow is absent. This makes it possible to accurately and unambiguously identify the onset and end of the periods with snow on the ground, in contrast to the conventional WMO synoptic observations, which are affected by occasional coding errors and often omit snow depth observations when snow is absent. However, some stations may not be representative of the general conditions prevailing in the surrounding area, owing to local effects. The HSDSD dataset starts in 1881 with a few stations and ends in 1995. A significant decrease in the number of available stations (from around 260 to around 190) occurred in January 1985 and in January 1994 (from 190 to less than 160). We considered the region extending from 35.5° to 73.5°N and from 30.5° to 180°W in order to cover the watershed of the major northern Eurasian rivers. During the time period considered, up to 263 stations provided daily snow depth observations, with approximately half of them covering the whole period without any missing data.
The second source of observations is the FSUHSS dataset (Krenke 2004). Data have been processed by the Institute of Geography, Russian Academy of Sciences (Moscow, Russia), and are accessible from the NSIDC data portal. It includes snow density, snow depth, and SWE observations. These observations were generally collected within one day, three times per month, along transects. Each snow depth observation is the average over 100 to 200 points around the station (less than 5 km), while SWE and density observations represent the average over 20 points. We used only the observations collected in open fields in order to ensure the consistency between both datasets. Density observations have been considered only when the snow depth was greater than 0.1 m, since density measurements are more difficult in the case of very shallow snowpacks. We verified that prescribing a threshold equal to 0.05 or 0.15 m does not significantly change the results.
c. Hydrometeorological datasets and reanalysis as sources of forcing data for snow models
Several global hydrometeorological datasets include meteorological variables at a 3-h resolution, which is sufficient to drive LSMs featuring detailed snow models. These variables are surface air temperature and humidity at a reference level (generally 2 m), surface wind velocity at a reference level (generally 10 m), surface incoming shortwave and longwave radiation, and precipitation. The simulation of the water budget over northern Eurasia is very sensitive to the forcing hydrometeorological dataset (Troy et al. 2011). Here, we used two different products: the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Re-Analysis (ERA-Interim) (Dee et al. 2011) and the Princeton University Global Forcing dataset (PGF) (Sheffield et al. 2006).
ERA-Interim meteorological variables have been extracted at 0.5° resolution. This is the latest generation of reanalysis produced by ECMWF and covers the time period from 1979 to the present (Dee et al. 2011). It is based on a 12-h four-dimensional variational data assimilation (4D-Var) system with a horizontal resolution of approximately 80 km and with 60 vertical levels. ERA-Interim incorporates several improvements with respect to the previous ECMWF reanalysis [40-yr ECMWF Re-Analysis (ERA-40)], which was based on a 3D-Var system. These enhancements include improved horizontal resolution (~80 km versus 125 km for ERA-40), the use of a variational bias correction to better handle satellite data biases (Dee and Uppala 2009), and the use of improved schemes to bias correct radiosondes and surface pressure observations. ERA-Interim also reflects advances in model physics, background error specification, data quality control, and radiative transfer modeling that make it possible to assimilate a greater number of satellite observations. An overview of the reanalysis can be found in Dee et al. (2011) and an evaluation of its performance is provided in Berrisford et al. (2011).
The PGF dataset is available at 1° resolution (Sheffield et al. 2006). It is based on the National Centers for Environmental Prediction (NCEP)–National Center for Atmospheric Research (NCAR) reanalysis. Sheffield et al. (2006) have performed corrections of the systematic biases in the 6-hourly NCEP–NCAR reanalysis via hybridization with global monthly gridded observations [National Aeronautics and Space Administration (NASA) Surface Radiation Budget (SRB) product for radiation, and Climate Research Unit (CRU) data for precipitation, temperature, and cloud cover]. In addition, the precipitation is disaggregated in both space and time at 1° resolution via statistical downscaling using the daily Global Precipitation Climatology Project data (GPCP). More details on this dataset can be found in Sheffield et al. (2006).
The major difference between PGF and ERA-Interim comes from the fact that the latter does not use any precipitation and radiation observations, either for the analysis of precipitation, or for the analysis of the other meteorological variables. ERA-Interim precipitation and radiation downward fluxes are those predicted by the meteorological model during the analysis process.
There are significant discrepancies between the existing datasets regarding climatological precipitation data in northern Eurasia (Troy et al. 2011). Therefore, building on former studies from Decharme and Douville (2006), we used Global Precipitation Climate Center (GPCC) gridded monthly precipitation (full data product V5) as an optional scaling factor for PGF and ERA-Interim precipitation.
We retained the period extending from 1 July 1979 to 30 June 1993, which combines the availability of ERA-Interim and PGF datasets together with a large number of stations providing snow depth records. To build a fair comparison framework between the simulations driven by the different datasets, we took into account the difference in elevation between each observation station and the nearest grid points. To achieve that, we applied the method developed by Cosgrove et al. (2003). This method is based on the standard altitudinal lapse rate of air temperature (−6.5 K km−1), which affects the altitudinal gradient of pressure accordingly. The specific humidity is changed under the assumption of a constant relative humidity and the longwave radiation is altered according to the changes in temperature and specific humidity, which affect the atmospheric emissivity. Wind speed, cloudiness, and total precipitation rates are assumed to remain unchanged, which is obviously a very simple assumption. In addition to the method developed by Cosgrove et al. (2003), the specific humidity has been adjusted when necessary to avoid a supersaturation with respect to ice in the case of subfreezing air temperature. For both datasets, we used the same snow–rain partitioning threshold (+1°C).
d. Snow model characteristics and setup
Simulations were performed with the Crocus snowpack model (Brun et al. 1989, 1992), now coupled with the ISBA LSM within the Surface Externalisée (SURFEX) interface (Boone and Etchevers 2001; Vionnet et al. 2012). SURFEX is the surface component of CNRM-CM (Voldoire et al. 2013). Crocus, now a snow scheme of ISBA among others, can be run either in a standalone mode, using a time series of prescribed meteorological observations as input data, or in a fully coupled mode with an interactive simulation of the atmosphere. The present study uses only the standalone mode. Figure 1 describes the main processes taken into account in the basic configuration. One important feature of the model is its ability to simulate snow metamorphism through a comprehensive set of semiempirical laws. The type and size of the crystals of each layer of the snowpack are prognostic variables, which control the surface albedo and the compaction rate of the different layers. Crocus has already been used in various research studies and operational applications. It has been widely evaluated, especially within the Snow Models Intercomparison Project (SnowMIP) (Etchevers et al. 2004).
The basic version includes a specific routine that detects the occurrence of blowing snow and empirically accounts for compaction and fragmentation of snow located in surface layers under wind-induced transport (Brun et al. 1997; Vionnet et al. 2012). However, mass losses due to blowing snow sublimation are not included. Blowing snow events occur when the local wind speed exceeds a threshold value determined by the state of surface snow (Schmidt 1982; Guyomarc’h and Merindol 1998). Such events may significantly affect the surface mass balance of snow-covered areas (Déry and Yau 2002; Yang et al. 2010). Indeed, wind transport creates heterogeneities in the snow cover through the formation of areas of erosion and accumulation. During airborne transport, snow particles also undergo sublimation when the atmosphere is unsaturated with respect to ice. This sublimation corresponds to a net mass loss for snow on the ground. To take into account blowing snow sublimation, a new process was added to Crocus as an option. Several parameterizations have been developed to estimate the sublimation rate in a column of blowing or drifting snow (Bintanja 1998; Déry and Yau 2001; Bowling et al. 2004). Gordon et al. (2006) combined some of them and suggested a simple parameterization for the total sublimation rate of blowing snow Qs, which has been implemented into Crocus. The quantity Qs depends on the near-surface meteorological conditions according to
where Ta is the air temperature (K), T0 a constant with a value of 273.15 K, U the wind speed, Ut the threshold wind speed for snow transport, ρa the air density, and Rhi the relative humidity with respect to ice. Here qsi denotes the saturation specific humidity (kg kg−1) at temperature Ta with respect to ice; γ, A, and B are dimensionless parameters; and Ut depends on the snow grain type at the surface simulated by Crocus (Guyomarc’h and Merindol 1998; Vionnet et al. 2012).
This parameterization makes the assumption that the change in the snowpack mass during wind-induced snow transport is controlled only by sublimation. It does not account for the possible local removal or accumulation of snow by wind, which is a frequent phenomenon depending on very finescale features of the wind field. Such a scale cannot be resolved in the context of the present study.
Since we aim to compare the snow model outputs with observations collected in open fields, the model is run assuming no interaction with vegetation and the snow fraction is set to 1, whenever snow is present on the ground. Simulations were performed using a version of ISBA featuring a multilayer explicit soil mass and heat transfer alternative to the original, simpler force–restore approach (Decharme et al. 2011). Heat flow proceeds along the thermal gradient, and soil ice is explicitly accounted for by the soil heat capacity, thermal conductivity, and soil hydrological parameterizations (Boone et al. 2000). Vertical soil water transfer is modeled using Richard’s equation, and the mixed form permits the modeling of a heterogeneous soil property profile. Soil freeze–thaw is explicitly modeled using a soil characteristic freezing curve, and soil ice is a prognostic variable. We prescribed a 10-m-deep soil profile, consisting of 11 numerical layers with increasing depths from 1 cm near the surface to 400 cm at maximum depth. Each layer uses a prescribed composition of ⅓ clay and ⅔ sand. We performed a first set of simulations starting on 1 July 1979 with a uniform soil temperature (5°C) and no snow on the ground. Then we used the soil temperature profile simulated on 30 June 1993 for each location as the initial state for a new set of simulations starting on 1 July 1979. This 14-yr spinup aims to initialize near-surface soil temperatures consistent with the local climate conditions of the different stations.
e. Statistical indices for evaluating the snow simulations
Various statistical indices were computed in order to evaluate the performance of the different simulations. More than 1 100 000 year-round daily snow depth observations have been used to derive the bias, correlation, and RMSE of the simulated daily snow depth. By taking into account only complete year-round series, amounting to approximately 2500 series, bias, correlation, and RMSe have been computed for the following annual variables: number of days per year with snow on the ground, duration, onset, and vanishing dates of the longest period with continuous snow on the ground (defined as snow depth higher than 1 cm). More than 100 000 observations of SWE and density have been used to calculate RMSE, bias, and correlation for these variables.
3. Results and discussion
Several simulations of the snowpack were performed, based on various combinations of the different options: PGF or ERA-Interim forcing, with or without monthly scaling of precipitation with GPCC data, and blowing snow sublimation on or off. Table 1 summarizes the performance of some of these configurations as a basis for discussing the following issues: the sensitivity to the driving dataset, the impact of the precipitation scaling, the impact of blowing snow sublimation, the reliability of ERA-Interim precipitation, the impact of the PGF precipitation disaggregation method, the evaluation of the simulated soil temperature, and the comparison with the Global Snow Monitoring For Climate Research (GlobSnow) product.
a. Results with the basic configuration of the snow model and the initial forcing datasets
Columns 1 and 2 of Table 1 show the performance of the basic snow model configuration driven with ERA-Interim and PGF, respectively. The simulations driven by ERA-Interim significantly outperform those driven by PGF. The performance is better for all statistical indices related to snow depth and the presence of snow on the ground. With ERA-Interim, the large positive snow depth bias is smaller than with PGF and all correlations are higher, revealing a better capturing of the intraseasonal, interannual, and spatial variability. Regarding SWE-related scores, RMSE and correlations are also higher with ERA-Interim, while the bias is slightly larger. This superior performance is a rather surprising result, since precipitation and radiation fluxes extracted from ERA-Interim are not based on observations but on meteorological forecasts, in contrast with the PGF dataset.
When computed over all dates, snow density scores are all significantly better with ERA-Interim than with PGF. This is also the case during almost the entire snow season, as shown in Fig. 2. With ERA-Interim, the model reproduces the observations and their standard deviations particularly well, with a difference in density generally smaller than 10%. The bias for the whole season is −6.6 kg m−3, which is about −3% of the mean observed density (222 kg m−3). The density simulated using the PGF forcing exhibits a significant negative bias until mid-April. This is partially due to less frequent strong winds in PGF than in ERA-Interim, probably induced by the lower spatial resolution of the PGF wind analysis. This results in a smaller compaction rate due to wind drift simulated by the model in this case (see section 3c). Toward the end of the season, especially in April, the density increases with time more in the model results than in the observations, leading to a good agreement between measured and simulated density data for both meteorological forcings. According to this comparison, it seems that Crocus overestimates the compaction rate in the presence of liquid water (i.e., during the melt season). This would explain why the density simulated using the ERA-Interim forcing exhibits a positive bias after mid-April. In contrast, the underestimation of the density simulated using PGF forcing would be partially mitigated during the melt season. However, the number of observations is much lower after mid-April than during the winter, which limits the reliability of the comparison between the simulated and the observed snow density during the melt season.
The spatial variability of density is also well reproduced with ERA-Interim, as illustrated in Fig. 3, which compares the observed and simulated midwinter snowpack density. The values correspond to the average density on 10 March from 1980 to 1993. An important point is the capacity of Crocus coupled with ISBA to accurately reproduce the density observed over eastern Siberia, which is around 150 kg m−3. This is a very low density since the snowpack in these regions starts to form in early October (6 months earlier). This performance stems directly from Crocus’ capacity to simulate the metamorphism of the snowpack layers into depth hoar under extreme temperature gradients. Such a realistic simulation of the snowpack density over a domain as large as northern Eurasia is a significant result of the present study. There are, however, a few stations, mainly in the northeastern part of Siberia, where the observed density is much larger than the simulation. This corresponds to stations that are strongly affected by accumulation and erosion induced by blowing snow.
b. Effects of precipitation scaling with GPCC data
Precipitation is obviously a major driver for snowpack simulations, which explains why the scaling of ERA-Interim precipitation with monthly GPCC precipitation has an important impact on performance (column 3 of Table 1). Most of the statistical indices are improved, except those related to density, which now exhibits a bias of −9.7 kg m−3 instead of −6.6 kg m−3. The improvement stems mainly from the drastic reduction in the snow depth and SWE bias, while maintaining high correlations. The simulation is almost unbiased in terms of SWE (+2.3 kg m−3; i.e., +3%) while the bias in snow depth is 0.008 m (+7%). The latter is explained by the overestimation of SWE and the underestimation of density (−9.7 kg m−3; i.e., −4% instead of −6.6 kg m−3). The annual snow duration and the dates of onset and melt out of the continuous snowpack show biases, which are less than half a week. In column 4 of Table 1, we notice a slight improvement of the simulations driven by PGF when precipitation is scaled with GPCC. All scores related to snow depth, snow cover duration, and SWE are improved, while those related to density are very similar. This confirms the superiority of GPCC compared to CRU over northern Eurasia in terms of monthly precipitation, as previously noted by Decharme and Douville (2006).
The significantly better performance of the snowpack simulations driven by ERA-Interim scaled with GPCC is not only true when averaged over the whole domain. Figure 4, based on the RMSE in the annual average snow depth, shows that this is also true for most individual stations. With the basic configuration of Crocus (no blowing snow sublimation), ERA-Interim with GPCC scaling gives the best performance over 154 stations, ERA-Interim without GPCC scaling over 69 stations, PGF without GPCC scaling over 20 stations, and PGF with GPCC scaling over 20 stations. The spatial distribution of the best-performing dataset exhibits no obvious geographical pattern.
The statistics in Table 1 summarize the performance over the whole set of stations. In general, the stations located in mountain or coastal areas exhibit poorer performances for both forcing data. For example, the bias in the annual snow duration is larger than 1 month for some of these stations. For this variable, as for the annual average snow depth discussed above, it is impossible to identify a regional pattern where PGF forcing outperforms ERA-Interim forcing, even when both datasets benefit from the GPCC precipitation scaling. The poorest performances in mountain and coastal regions probably stem from the higher spatial variability of the meteorological conditions, which is not captured by either meteorological dataset.
c. Effects and representation of blowing snow events on snowpack sublimation
Column 5 of Table 1 shows the results of the snowpack simulation driven with ERA-Interim when blowing snow sublimation is taken into account. When compared with the basic configuration (column 1), the performance is greatly improved for snow depth and SWE, both in terms of bias and RMSE. The annual snow duration and the dates of onset and melt out of the continuous snowpack also exhibit improved biases (2 days or less in absolute value). The bias in snow density is worse (−10.3 kg m−3) but is only −5% of the average observation. The simulation with blowing snow sublimation and the simulation without this parameterization but scaled with GPCC (column 3 of Table 1) perform similarly.
The drastic reduction of the SWE and snow depth biases is in agreement with previous studies. Pomeroy and Gray (1995) estimated that losses due to blowing snow sublimation may reach 15%–41% of annual snowfall over the Canadian Prairie. In arctic Alaska, losses were found to range between 9% and 22% of annual snowfall by Liston and Sturm (1998). At larger scales, Yang et al. (2010) studied the effects of blowing snow transport and sublimation on the winter snow season mass balance of the Northern Hemisphere, using a mesoscale atmospheric model coupled with a blowing snow model that includes sublimation. They found a weak contribution of snow transport at large scales, but identified local maxima of blowing snow sublimation in the northwest of Mongolia, central Kazakhstan, and northern Siberia, where open fields dominate.
Taking into account the sublimation associated with blowing snow events improves not only the global statistics, but also the distribution of the statistics for the individual stations as well. Figure 5 presents the distribution of the bias in the date of the end of the continuous snow cover. Introducing the blowing snow sublimation leads to a much better symmetry between the number of stations showing negative and positive biases, and more than half of the stations exhibit an absolute bias less than 6 days. In spite of this obvious improvement, there are still seven stations where the bias exceeds 1 month. This mainly corresponds to stations where the local winds are much stronger than ERA-Interim winds, leading to frequent and strong blowing snow events. The corresponding snowpacks are very shallow and vanish more rapidly than in the simulations.
To understand the effect of the parameterization of blowing snow sublimation, we looked into details of the individual simulations. Each station features a distinctive behavior, but some patterns are common to most of them. As an illustration, Fig. 6 compares the daily snow depth observed at Ivdel station (60.7°N, 60.4°E) with the simulation performed with and without the blowing snow sublimation, respectively. The figure also shows the average ERA-Interim daily wind speed and the corresponding observation. It illustrates the principal impacts of the parameterization:
The introduction of blowing snow sublimation limits the snow depth overestimation. Even if the observed and the analyzed ERA-Interim wind speeds frequently differ, they are rather well correlated and the parameterization significantly improves the simulations for many cold seasons, especially in 1980/81, 1981/82, 1983/84, 1985/86, 1986/87, 1988/89, and 1990/91.
Details show that the benefit of the parameterization is due to very few blowing snow events that induce a significant decrease both in the simulated and in the observed snow depth. The vertical dashed lines in Fig. 6 highlight the correspondence between relatively strong winds and decreases in simulated and observed snow depth. The differences induced by these events persist for weeks and even months. This explains a large part of the bias found in the version of Crocus without blowing snow sublimation.
Some winters show only small differences between both simulations (for example winter 1991/92, which exhibits very few windy situations). This illustrates the interannual variability of blowing snow events, which depends on the occurrence of high wind speeds together with transportable/erodible layers at the top of the snowpack.
The simulation with the parameterization is poorer during 1984/85, with a slight underestimation of snow depth. That is mainly due to a unique blowing snow event around mid-December, which is simulated, but not observed.
d. Reliability of ERA-Interim snow precipitation over northern Eurasia
Among the different configurations, two simulations exhibit similar scores for most of the statistical indices (columns 3 and 5 of Table 1). Both simulations are driven by ERA-Interim, but the small bias is obtained either from the precipitation scaling with GPCC or from accounting for blowing snow sublimation. The simulation including the precipitation scaling together with the blowing snow sublimation leads to a significant negative bias (column 6 of Table 1). There is an obvious redundancy between scaling ERA-Interim precipitation with GPCC and taking into account blowing snow sublimation.
Indeed, the total amount of ERA-Interim precipitation over December, January, and February is 18% higher than GPCC precipitation when averaged over the considered stations and over the simulation period. But GPCC uses uncorrected precipitation gauges and is therefore recognized to underestimate snowfall. In our domain of interest, this underestimation can reach up to 25% (see GPCC correction over the winters 2007–10 at http://kunden.dwd.de/GPCC/Visualizer). It is linked to gauge undercatch during windy conditions (Sevruk 1987; Legates and Willmott 1990; Goodison et al. 1997). In addition, the work by Troy et al. (2011) on the terrestrial water budget in northern Eurasia did not reveal any significant bias in ERA-Interim winter precipitation. Therefore, there is no strong suspicion of a possible overestimation of ERA-Interim snowfall.
To consolidate this assertion, we have evaluated the surface snow sublimation simulated by ISBA/Crocus during the winter. By surface snow sublimation, we refer solely to the turbulent mass exchanges between the surface of the snowpack and the atmosphere, excluding the sublimation associated with blowing snow events. It depends on snow surface temperature, air temperature, air humidity, and wind velocity. Figure 7 shows the 2D field of sublimation/condensation amount simulated by Crocus during winter [December–February (DJF)] over open fields, averaged from 1979 to 1993. These amounts have been compared with estimations from former studies that have addressed this issue. The comparison makes sense only over areas where vegetation is sparse since our simulations are valid only over open fields. Yang et al. (2010) have estimated winter surface sublimation over DJF 2006/07 with a coupled LSM–atmospheric model. Over the nonforested areas of northern Siberia, they calculate condensation rather than sublimation with values varying from 0 to +5 kg m−2, which is in agreement with our estimations. With the LSM VIC, Troy et al. (2011) have estimated evapotranspiration rates lower than 1 or 2 kg m−2 per month during DJF over the major Siberian watersheds. However these watersheds include forested regions and the results over the nonforested parts are not distinguished, which limits the comparison with our simulations. Nevertheless, their estimations are not significantly higher than the maximum sublimation rates that we have simulated over open fields in the western part of our domain (typically from 0 to −5 kg m−2 over DJF). These considerations make us confident in our simulations, which show that the DJF surface sublimation contribution to the snow cover mass budget is very low over northern Eurasia (condensation or sublimation from +5 to −5 kg m−2 from northeast to southwest). Consequently, a possible underestimation of the simulated sublimation could not explain the positive bias in snow depth and SWE that the basic configuration of Crocus shows when driven by ERA-Interim without accounting for blowing snow sublimation.
All the facts invoked above work together in favor of a very good reliability of ERA-Interim snow precipitation. The redundancy between scaling ERA-Interim precipitation with GPCC and accounting for blowing snow sublimation has a physical explanation. The meteorological conditions that lead to an undercatch of snowfall by precipitation gauges also correspond to periods inducing blowing snow, since fresh snow layers exhibit a low density and thus snow can be transported even under moderate wind speeds. The frequent simultaneity of both processes during a snowfall explains the fact that scaling ERA-Interim with GPCC is a subterfuge for taking into account the blowing snow sublimation. This also explains why scaling with GPCC and accounting for blowing snow sublimation lead to a significant negative bias (column 6 in Table 1).
e. Identification of PGF temperature/precipitation inconsistency
The snow cover simulations with the basic configuration of Crocus and with forcing data from ERA-Interim or PGF (columns 1 and 2 in Table 1) are all characterized by a large positive bias in snow depth and SWE (20 kg m−2; i.e., 29%). However, ERA-Interim simulations show significantly higher correlations, and scaling ERA-Interim precipitation with GPCC reduces this bias considerably (column 3), which is not the case with PGF (column 4). That means that the larger bias with PGF is due to meteorological variables other than the amount of precipitation. To identify whether the best performance of ERA-Interim simulations is due to one specific meteorological driving parameter, we performed additional simulations with forcing data obtained from a hybridization of PGF and ERA-Interim temperature/humidity and precipitation. ERA-Interim air temperature and humidity give better results than their PGF counterparts. PGF precipitation, whether scaled with GPCC or not, improves the snow depth bias but decreases the correlation. To attribute a cause to this behavior, we looked at individual stations in detail. While there are more than 2500 samples (1 station, 1 year) and contradictory examples coexist, we consider the simulations for Pirovskoe (57.6°N, 92.3°E) in Fig. 8 to be a representative illustration of a behavior commonly noticed for many stations. The snow depth simulated with PGF is overestimated during several weeks and even months during the cold seasons 1979/80, 1981/82, and 1982/83, mainly as a result of a few early snowfalls which occurred around mid- to late October. Figure 9 shows the chronology of the ERA-Interim, PGF, and observed daily precipitation and the simulated snow depth during October 1982. Indeed, in the observations, as in ERA-Interim, the most intense precipitation actually fell during the relatively warm period 4–7 October. In contrast, the significant precipitation in PGF during the relatively cold period 9–17 October is not observed. ERA-Interim undoubtedly exhibits a much better precipitation chronology, resulting in a better snow–rain partitioning in the simulations.
This example, representative of many other stations not shown here, clearly demonstrates that the ERA-Interim chronology of precipitation during the cold season is much better than its PGF counterpart. In PGF, monthly snowfall events are generally distributed over very few days, which explains the frequent large peaks in the snow depth simulation that are not observed (see Fig. 8). The better quality of the ERA-Interim snowfall chronology is explained by the fact that cold season precipitation over northern Eurasia is mainly caused by synoptic-scale systems, with mostly stratiform precipitation. These systems are well captured in ERA-Interim. The poorer simulation results obtained with PGF are an unexpected consequence of the monthly precipitation temporal disaggregation method described in Sheffield et al. (2006). This method distributes the monthly precipitation over a limited number of days within the month in order to better match the climatology. This is problematic for snowpack simulations, as the temperature of these days may differ significantly from the temperature of the days during which precipitation actually occurred, changing the snow–rain partitioning. The disaggregation method in PGF has been recognized as a major advance for many hydrological studies but it has a negative impact on snowpack simulations, particularly at the onset of the snow season when temperatures are still often close to the melting point. This period of the year generally exhibits heavier precipitation than midwinter, which explains its significant impact on snow depth and SWE during most of the snow season.
To the best of our knowledge, the negative impact of the temporal precipitation disaggregation in PGF on the snow depth simulations has never been noticed before. This is probably because it occurs in regions with low precipitation and because the total amount of available water is unchanged when we add rainfall to snowmelt runoff over the course of the cold season. However, it does have a significant impact on the simulated snowpack characteristics and limits the reliability of snow cover simulations driven by PGF.
f. Evaluation of near-surface soil temperature
We have attributed the good performance of the density simulation to Crocus’ ability to simulate the formation of depth hoar under high temperature gradients. This type of snow crystal inhibits the compaction of the snowpack and maintains a low snow density. Simulating large temperature gradients in the snowpack requires the simultaneous presence of cold snow surface temperatures and relatively warm soil surface temperatures. We have verified that the simulated snowpack consisted primarily of depth hoar snow layers. However, this study uses the recent coupling of Crocus with the multilayer soil model ISBA soil model, which includes freezing-soil processes. An evaluation of the simulated near-surface soil temperatures is necessary, especially in regards to the relatively short spinup (14 years) used in our configurations to initialize the soil component of the model. We have compared the simulated temperature at 0.20-m depth with the monthly observations available in the Russian Historical Soil Temperature Data dataset (Zhang et al. 2001). We used only those stations that are collocated with the HSDSD or FSUHSS stations and that exhibit a difference in elevation less than 100 m from the corresponding meteorological forcing grid points. During the considered period, the evaluation dataset amounts to 96 stations, for a total of 11 760 monthly observations. Figure 10 presents the results for the simulation driven by ERA-Interim and accounting for blowing snow sublimation. With a bias equal to −0.5 K, a correlation equal to 0.675 and an RMSe equal to 2.8 K, the simulated soil temperature at 0.20-m depth is quite realistic. Since the simulated near-surface soil temperature is very sensitive to the snowpack characteristics (Gouttevin et al. 2012), this performance is an additional indirect evaluation of the snowpack density and depth simulated with Crocus and an evaluation of the capacity of the multilayer soil model ISBA to simulate the soil freezing and thawing.
g. Comparison with GlobSnow products
We computed the performance of SWE simulations over all stations where survey transects are available (1089 stations during the considered period) and compared it with the performance of GlobSnow retrievals (Luojus et al. 2011). The domain and the period are not exactly the same (35°–85°N for GlobSnow instead of 35°–73.5°N and 1979–2000 instead of 1979–93). Table 2 compares the bias and RMSE on SWE from GlobSnow and Crocus against observations. As in Luojus et al. (2011), a distinction is made between the statistics derived from all observations and those derived from the subsample limited to observations with a SWE lower than 150 kg m−2. Results in Table 2 show similar performances for both products. This is very encouraging for the snowpack simulations since they do not incorporate any snow depth or SWE observations, in stark contrast to GlobSnow products, which assimilate the snow depths recorded at synoptic stations.
Using ISBA/Crocus, a detailed snowpack model coupled with a soil model, and forcing meteorological data from ERA-Interim and PGF, we have performed snow simulations and compared them with daily snow depths observed over 14 years across more than 250 stations in northern Eurasia. The best simulation is obtained with ERA-Interim. Without any local calibration or integration of snow depth or SWE observations, the bias in snow depth is 0.006 m (i.e., about +5% of the average value) and the bias in snow duration is only −1.2 days. The bias in the dates of the onset and end of continuous snow cover is very small as well (−1.2 days and 2.0 days, respectively). Correlations are also quite high, revealing a reasonable capturing of the spatial and interannual variability.
Regarding local SWE observations, the performance of the snowpack simulations is very similar to the accuracy of gridded snow products from satellites which assimilate local snow depth observations. The spatial and temporal variability of snow density is also well simulated, thanks to the capacity of the snow model to simulate the formation of depth hoar and consequently the considerably lower compaction rates. This performance stems from the capacity of the coupling between Crocus and the multilayer soil model ISBA to simulate the snowpack thermal resistance and the soil freezing and thawing reasonably well.
We have shown the importance of blowing snow sublimation and the parameterization from Gordon et al. (2006) proved to be very efficient. We have shown that there is a redundancy between accounting for this process and scaling ERA-Interim precipitation with monthly precipitation datasets based on uncorrected precipitation gauges, like GPCC. This is due to the fact that the windy conditions leading to an undercatch of snowfall by precipitation gauges induce blowing snow events as well. Therefore, we recommend scaling ERA-Interim precipitation with GPCC only for those models that do not simulate the blowing snow sublimation.
The methodology of the present study is also an original way to compare the reliability of different meteorological datasets. It revealed that the PGF dataset suffers from an inconsistency between the chronology of precipitation and temperature, limiting its ability to simulate a realistic onset of snow cover. That affects the simulated snow depths during a significant part of the winter. We recommend that future studies on snow modeling over this region take this issue into account, especially if they focus on the onset of snow cover. An updated version of PGF correcting this inconsistency would be very useful. We have shown that ERA-Interim snowfall is very reliable. The intrinsic consistency between all its meteorological variables and its general performance favor a realistic simulation of the energy and mass budget of the snow cover over open fields. This result is consistent with previous works on ERA-Interim radiation (Troy and Wood 2009) and makes ERA-Interim a very attractive dataset for snow modeling.
To the best of our knowledge, the present study is the first to evaluate, in detail, snowpack simulations over open fields against more than 1 100 000 local daily snow depth observations, and over more than 100 000 SWE and density observations that have been measured over open fields as well. This is a major step in snow modeling: it proves that meteorological reanalysis and LSMs, including detailed snow models, have now reached a performance level sufficient to reconstruct key features of open field snowpacks over very large regions without any local calibration. Should our study be extended to the whole of North America and Eurasia, snow density simulations could be very useful for real-time snow products that assimilate open field snow depth observations from synoptic stations. Indeed, it would provide a reliable source of information for transforming snow depth observations into SWE estimations according to the local climate conditions and complement recently developed approaches (Sturm et al. 2010). It could also improve the estimation of SWE from satellite-based snow retrievals algorithms (Boone et al. 2006). ISBA/Crocus simulations also include a comprehensive set of internal snowpack characteristics like the temperature, density, grain size, and shape of the different layers it is made of. Such variables could be very useful for the improvement of remote sensing algorithms.
Finally, the quality of the simulations over open fields is very encouraging for future implementations of detailed snow models into GCMs. This is also the case for the development of real-time snow cover monitoring systems based on operational meteorological analyses and detailed snow models, with or without the assimilation of snow depth observations. However, the performance of such applications requires a reliable simulation of snowpack properties under the forest canopy. The present study is limited to snowpack simulations over open fields, which is easier than in the forests. Crocus has no physically based option to take into account the effects of vegetation at the present time, but it will be implemented in the near future. SWE and density observations are available from some forest stations in the FSUHSS dataset, offering the possibility of extending the present study from open fields to the forest.
This work was funded by ANR, under CLASSIQUE project number ANR-2010-012-03, and by the European Commission’s 7th Framework Programme, under Grant Agreement 226520, COMBINE project. The authors are thankful to the Institute of Geography, Russian Academy of Sciences (Moscow, Russia), the State Hydrometeorological Service (Obninsk, Russia), and NSIDC (Boulder, CO), which collected, processed, or delivered the snow observation data used in the present study. They are also thankful to the numerous observers from the former Soviet Union who performed these precious observations, often under very harsh meteorological conditions. The authors thank Paul Poli, from ECMWF, who provided detailed information on ERA-Interim assimilation processes.