Evaluating global snow water equivalent products for testing land surface models

This paper compares three global snow water equivalent (SWE) products, SSM/I (NSIDC), AMSR-E (NSIDC) and Globsnow (v1.0, ESA) to each other, snow covered area (SCA), ground measures of snow depth and meteorological data in an attempt to determine which might be most suitable for testing and developing land surface models. Particular attention is paid to which gives the most accurate peak accumulation, seasonal SWE changes and first and last dates of snow cover.    SSM/I and AMSR-E are pure earth observation (EO) derived products whilst Globsnow is a combination of EO and ground data. The results suggest that the pure EO products can saturate in deeper snow (SWE > 80–150 mm), can show spurious features during melt and can overestimate SWE due to strong thermal gradients and erroneous forest cover correction factors. Along with the comparison to ground data (only a single point) this suggests that Globsnow is the more accurate product for determining peak accumulation and seasonal SWE cycle.    The snow start and end dates of the three SWE products were compared to an optically derived SCA (MOD10C1, taken as truth) and found to give large errors of snow start date (root mean square error of 20 + days, though SSM/I was correct on average). The snow end dates had lower errors (a bias of 1–6 days) although the spread was still on the order of three weeks.    During the investigation, occasional abrupt changes in Globsnow were observed (in the v1.0 and v1.2 daily and v1.2 weekly products). These only occurred in around 1% of cases examined and seem to be spurious. Care should be taken to corrector avoid these jumps if using Globsnow to validate land surface models or in an assimilation scheme.

section_content: Snow is the largest transient feature of the land surface (Yang et al. 2001). It has an effect on atmospheric circulation through changes to the surface albedo, thermal conductivity, heat capacity and aerodynamic roughness (Gong et al. 2004), as well as controlling the availability of water in some ecosystems and to a sixth of the world's population (Clifford 2010). Therefore it is vital that snow is properly represented in land surface models if we are to understand and make predictions of weather, climate, the carbon cycle, flooding and drought. 

To test snow processes in land surface models (LSMs) global measures of snow are needed; in particular the presence or absence of snow (snow covered area, SCA) and the amount of snow (snow water equivalent, SWE). Satellite earth observation (EO) is the only way to collect globally consistent data regularly, especially as many snow-affected areas are sparsely populated (Davenport et al. 2012). EO estimates require inversion algorithms to relate raw signals recorded at the satellite to physical properties of the land surface and these inverted estimates can contain errors and biases. 

Whilst EO SCA is considered reliable (Hall and Riggs 2007, Brown et al. 2010), SWE has proven more problematic, especially during melt (see section 2. 1). Some studies suggest only trusting EO SWE measurements up to peak accumulation, relying on SCA alone after that (Déry et al. 2005, Kuchment et al. 2010). However, others have found consistent differences between EO peak SWE estimates and other data sources, although which was correct was not known (Clifford 2010). 

This paper aims to test three global, EO-derived SWE products against EO-derived SCA, meteorological data and some limited ground measurements to attempt to identify any issues or biases that would impair their suitability for validating land surface models.
section_name: Introduction

section_content: The three SWE products tested were NASA/JAXA's AMSR-E/Aqua Daily L3 Global Snow Water Equivalent EASE-Grids (AE DySno) (Tedesco et al. 2004 ), NSIDC's Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover (NSIDC-0321) (Brodzik et al. 2007 ) and ESA's Globsnow SWE v1. 0 (Takala et al. 2011) ; hereafter referred to as "AMSR-E", "SSM/I" and "Globsnow". All three cover at least the northern hemisphere at 25 km spatial resolution although the standard Globsnow product masks out mountainous areas (see section 2. 4 for details), a version with estimates over mountains is available, though this is unvalidated. Unfortunately mountainous regions are important for water supply in some areas (Clifford 2010). AMSR-E and Globsnow provide daily estimates (with occasional gaps towards the mid latitudes), whilst SSM/I provides a measure every eight days (an average, though how the multiple measurements and platforms are combined is not obvious from the product description). The theory of measurement and inversion algorithms are described below. Since this work, Globsnow v1. 2 and v1. 3 have been released and section 3. 7. 1 briefly examines v1. 2. v1. 3 was released during the paper's revision and is the same as v1. 2 except over North America.
section_name: Data used

section_content: SSM/I and AMSR-E are passive microwave radiometers. All objects emit microwaves, although the soil emits a stronger signal than the snow-pack and these microwaves are scattered and attenuated by the snow as they travel upwards. The attenuation will depend on the microwave wavelength and the snowpack properties, in particular the amount of snow (SWE), grain size, snow density, presence of ice lenses and the amount of liquid water (Ulaby et al. 1982). The shorter the wavelength the greater the scattering for a given set of snow-pack properties, thus the difference between the signal at two wavelengths should be related to the amount of snow. This has been theoretically derived and tested over certain areas (described in section 3. 6) to give the following equation (Foster et al. 1980) ;

Where SD is snow depth in cm, T b 18H is the brightness temperature at 18GHz, horizontal polarisation, T b 36H is the brightness temperature at 36 GHz, horizontal polarisation, and a is a constant (generally set at 1. 59). This can be converted to SWE by multiplying by snow density (assumed to be 300 kgm -2 ). This is known as the "Chang algorithm" (Chang et al. 1987) and forms the basis of the AMSR-E and SSM/I inversions. The wavelengths will vary for different instruments and correction factors can be applied to account for this (Brodzik et al. 2005). 

Although it has been validated to some extent, this approach is known to suffer from a number of issues. It has been shown that snow grain size and SWE dominate the scattering (Chang et al. 1982). Snow density has an effect, although this is much smaller. The larger the snow grains the greater the scattering for a given SWE and so a, in equation 1, is tuned for a particular snow grain size. The value of 1. 59 was found by assuming a grain size of 0. 35 mm (Chang et al. 1982), but as snow grains change through the winter an assumption of constant grain size will cause errors (Kelly et al. 2003). Liquid water affects the snow-pack's dielectric constant, increasing attenuation, potentially preventing the retrieval of SWE (Stiles and Ulaby 1980). Both AMSR-E and SSM/I take measurements at night when it is hoped that the snow has refrozen (AMSR-E crosses the equator at 1:30 am and SSM/I at 6 am, within 20 minutes depending on which platform (NSIDC 2012), Globsnow uses these same sensors, see section 2. 4); however, any melt and refreeze is likely to alter the grain size and could lead to errors. 

The Chang algorithm has also been shown to saturate in deeper snow covers, SWE >120 mm (Chang et al. 1987). This is due to the higher frequency microwaves no longer decreasing with increasing SWE (Derksen et al. 2010). Therefore accuracies are reduced in deep snow. 

As well as the snow, any vegetation will also scatter the microwaves, affecting the estimate. A forest correction factor, c, has been proposed (Foster et al. 1997). 

Where f is the fractional forest cover within a pixel. Whilst this makes sense theoretically, there is some doubt as to its accuracy (Clifford 2010) and it may lead to inaccurate estimates over forested areas. 

It has been suggested that topography can introduce errors (Brodzik et al. 2007), both through changes in the atmospheric thickness and because AMSR-E and SSM/I scan off nadir. Finally, shallow snow packs cause little scattering and can be missed by microwave sensors; a minimum detectable SWE of 15 mm has therefore been suggested (Solberg et al. 2010 ).
section_name: Microwave measurement of snow

section_content: SSM/I is a passive microwave radiometer onboard the US Defense Meteorological Satellite Program polar orbiting platforms, a number of which have been in orbit from 1988 to the present day (F8, F11, F13 and F17 (NSIDC 2012) ). The SWE product, as provided by the National Snow and Ice Data Centre (NSIDC), extends from 2000 to 2008 and is available for download (Brodzik et al. 2007). This product uses the Chang algorithm, as given in equation 1, multiplied by a forest correction of the form given in equation 2. The forest correction factor has a maximum value of 2 for 50% and above forest cover, the latter being derived from the BU-MODIS land cover map (Knowles 2004). See section 2. 8 for a description of the BU-MODIS dataset.
section_name: SSM/I

section_content: AMSR-E is a passive microwave radiometer mounted on NASA's Aqua satellite. It has been collecting data since launch in 2002 but as of the 4 th October 2011 a fault has stopped measurements (NSIDC 2011). The SWE product is available from NSIDC (Tedesco et al. 2004 ). 

The AMSR-E algorithm is described in detail in Kelly et al. (2003) and Kelly (2009). This is based upon the Chang algorithm, but the grain size can change with time. Fresh snow grains start at 0. 2 mm and grow at a rate determined by surface temperature (derived from AMSR-E polarisation differences) until they reach 1 mm. This is then used to tune a slightly more complex version of equation 1 (still depending on the difference between the brightness at different frequencies, though with the addition of 10 GHz to help with forest attenuation). This is quite a simple grain growth model, but it was hoped it would lead to more accurate SWE estimates. A forest correction is also applied using the product of the MOD12Q1 and MOD44B products, both derived from MODIS (Barnes et al. 1998). See section 2. 8 for a description of MOD12Q1. MOD44B is derived from MODIS and gives the proportion of forested, non-forest vegetated and bare ground within a 250 m pixel (Hansen et al. 2006). 

This method still depends on the same basic assumption as the Chang algorithm and may suffer from similar issues.
section_name: AMSR-E

section_content: In August 2010 ESA released the Globsnow v1. 0 SWE and SCA products (Takala et al. 2011 ) and these are available online (http://www. globsnow. info/). Globsnow is a fusion of EO-derived estimates with ground measurements and covers the period from 1978 to the present day. The most modern passive microwave radiometer available is used (SMMR 1978 -1988, SSM/I 1988 -2002, AMSR-E 2002 -2011 and SSM/I since October 2011). Ground measurements of snow depth are taken from the ECMWF WMO weather stations (Takala et al. 2011). These are point measurements (or snow courses on the order of a few metres). 

Using the method described in Pulliainen (2006), the HUT (Helsinki University of Technology) microwave snow emission model (Pulliainen et al. 1999 ) is used to relate the brightness temperature values of the satellites to the snow properties and snow depth. HUT is a radiative transfer model for calculating the microwave emission of a snowpack (assumed homogeneous) for given SWE, grain size and snow density. By inverting HUT for the observed brightness temperature (from the satellite), snow depth (average of the six closest ground stations) and fixed snow density (240 kgm -2 ) an effective microwave grain size is derived. The standard deviation for the six ground stations gives a measure of the uncertainty. This uncertainty is combined with a HUT derived sensitivity to changing SWE so that uncertainty increases as the sensitivity decreases (as SWE approaches the saturation point). 

These grain size values are interpolated by kriging to produce a global map. This is used to invert HUT with a fixed grain size and observed brightness temperature to give an EO estimate of snow depth; converted to SWE by assuming a density of 240kgm -3. The ground depth measurements are converted to SWE and interpolated to give a global ground map. This is assigned a variance of 250 mm 2, based on comparison of point depth measurements to long snow courses (Takala et al. 2011). 

These two estimates are combined by minimising the difference between the observed brightness temperature and snow depth for a given snow depth, weighted by their uncertainties, to give a global SWE map along with an uncertainty estimate. The grain size is allowed to deviate from the initial estimate during this minimisation and the EO uncertainty is adjusted accordingly. Additional dry snow start (Hall et al. 2002) and snowmelt start (Takala et al. 2009 ) algorithms are employed. Before the first dry snow and after the onset of melt the EO weight is set to zero. 

In theory this should provide more accurate estimates of SWE than either of the Chang algorithm products (AMSR-E and SSM/I), as the grain size is derived from ground data and forward snow emission model simulations with direct measurements of depths are used. An accuracy of ±40 mm SWE (RMSE) was the aim, though this may increase for SWE above 150 mm due to challenges and uncertainties when modelling deep snow packs (Takala et al. 2011). 

Ground stations tend to be in valleys rather than on mountain tops and so may not be representative of the conditions within a topographically complex pixel. For this reason (and the issues with microwave data over complex topography) all 25 km pixels with an elevation standard deviation greater than 200 m are marked as mountainous and no SWE estimate is provided (Takala et al. 2011) ; the ETOPO5 terrain map is used to derive elevation range. In September 2011 an unmasked version was released, though this includes no extra processing and estimates over mountainous regions may not be accurate (Luojus, pers comms). Version 1. 0 was released in March 2011 and was used for this study. Version 1. 2 was released in October 2011 and version 1. 3 in August 2012; this uses the same theory as above but with slight differences in the assimilation and use of ground data. 

As well as SWE, the grain size maps derived during Globsnow's creation are recorded (daily, 25 km resolution) and are available on request from the Globsnow team. The locations of ground stations contributing to Globsnow each day are also recorded and available, although the weightings of these and the EO data are not.
section_name: Globsnow

section_content: Two snow covered area products were used in this study, MOD10C1, full name "MODIS/Terra Snow Cover Daily L3 Global 0. 05Deg CMG (Hall et al. 2006 ) and the MODIS-SSM/I blended product, full name "NSIDC-0321 Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover" (Brodzik et al. 2007), hereafter referred to as MOD10 and MODIS-SSM/I. MOD10 is derived from the MODIS passive optical sensor (2000-present) (Barnes et al. 1998). It uses the "Normalised Snow Difference Index" (NDSI), which is the ratio of the difference in reflectance in the infra-red and the visible to the sum of the two, to estimate the fractional snow cover within a 500 m pixel once a day (Salomonson and Appel 2004). This only relies upon snow changing the spectral properties of the surface. Therefore, it should give reasonably accurate estimates as long as clouds do not obscure the ground, irrespective of snow-pack structure and water content. Hall and Riggs (2007) report an accuracy compared to ground observations of 93% on average, though over evergreen forests this decreased to 80% due to shadowing effects (Painter et al. 2009). 

Cloud cover can be a persistent problem in some areas and optical measurements are not possible during the Arctic winter night. To overcome this the passive optical product can be supplemented with passive microwave estimates, providing a continuous SCA estimate (Liang et al. 2008). A blended product uses MODIS when available and a method similar to that of Chang et al. (1987) with SSM/I when it is not. The microwave sections will suffer from the same issues with liquid water, snow-pack structure and shallow snow layers as the SWE estimates described above and so may not be as reliable as MOD10. In addition the MODIS-SSM/I product has the added advantage of being 25 km resolution, rendering global analyses more practical.
section_name: Snow cover

section_content: The SWE estimates described above are sensitive to snow-pack structure that in turn is controlled by air temperature. In addition, it may be possible to use cumulative snowfall as an approximation for SWE whilst the air temperature stays below 273K. Cumulative snowfall will not exactly match SWE due to sublimation and mid-winter melt (redistribution should not be an issue at the scales investigated here), but it will be related and so provide a guide. Global estimates of both can be provided by re-analysis data (a combination of general circulation model, GCM, runs with assimilated ground measurements to reduce any bias). For this study the dataset of Sheffield et al. (2006), hereafter referred to as "Princeton", was chosen because it overlaps with all of the above snow products. 

Princeton provides (amongst other variables not relevant here) estimates of air temperature and precipitation (total rain and snow) every 3 hours from 1948 to 2008 at 1 o spatial resolution. This is a much coarser spatial resolution than those of the snow products; all data were therefore aggregated to this by simple averaging before comparison (other than for a few tests performed at the native resolution). As the EO data should be representative of the properties across each pixel this should not introduce any error. It is also the resolution at which LSM validation would take place. 

The air temperature should be reliable but there are great problems measuring precipitation in winter. Snow is easily blown, missing precipitation gauges, and can also clog them (Groisman et al. 1991). Small snowfalls can be difficult to detect and mixed snow and rainfall can complicate attempts to correct for wind (Yang et al. 1998). Therefore Princeton may underestimate the snowfall. At more intensively managed stations this should not be as much of an issue, but it is always a possibility. Therefore, we should not rely on absolute snowfall values but use it as a rough guide. 

The total precipitation provided by Princeton can be partitioned into rainfall and snowfall using a temperature threshold, typically between 273K and 275K (Best 2005). It was set to 274K for this study.
section_name: Meteorology

section_content: Ground measurements of SWE should be much more accurate than EO and re-analysis estimates; however, they only cover limited areas and time periods and may not be representative at the scales of the products described above. 

An automatic weather station was maintained near Kevo, Finland between November 2007 and April 2009 (69 o 45'N, 27 o 1'E) as part of the ABACUS project. This included a sonic snow depth measurement using a Campbell SR50 Ultrasonic distance sensor (Campbell Scientific, Logan, Utah, USA) every thirty minutes. This was converted to SWE assuming a density of 300 kgm -2. This density was taken from the literature (Chang et al. 1987) and was not tuned to match the datasets in any way. 

The gauge works by recording the time of flight for a click to travel from the sensor to the snow and back. This click will bounce back from the first object, whether that is the ground, a blown snowflake, tuft of grass or a passing hare. Spurious results from such encounters were obvious in the data and were cleaned by smoothing with a five-day median filter followed by convolution with a three-day Gaussian. 

This was a temporary weather station and is not used in the creation of Globsnow. The nearest ground station used was approximately 50 km away.
section_name: Ground data

section_content: Snow processes and the measurement of snow are also affected by topography and forest cover. Global measures of topography were taken from NASA/JAXA's ASTER Global Digital Elevation Model version 1 (G-DEM), which is available online (NASA 2009). This gives global elevation at 30 m horizontal resolution, however it was coarsened to 150 m with a median filter to remove artefacts as recommended by Hayakawa et al. (2008). 

Forest cover was calculated from the BU-MODIS land cover map (Knowles 2004 ). The BU-MODIS product is a flat binary version of NASA's MOD12Q1. That is a global land cover map derived from MODIS data collected between the 1 st January 2001 and 31 st December 2001. This uses optical information and the changes over a year to class each 1 km pixel as one of the seventeen International Geosphere Biosphere Programme land cover classes (Strahler et al. 1999). These seventeen classes were used to determine the forested fraction of each 1 o pixel.
section_name: Other datasets

section_content: This section describes the experiments performed and the areas they were performed over. As the methods rely on the results of earlier sections, the method, results and discussion for each experiment are grouped.
section_name: Experiments and Results

section_content: A subset of 1,381 1 o pixels was used; chosen to be representative of the global snow-affected areas in terms of latitude, longitude, elevation, topographic variation and peak annual SWE (hence the unusual number). The number of snow-affected land pixels (taken as pixels marked as snow covered on the 25 th March 2005, an arbitrary date close to peak accumulation) was counted and stratified random sampling performed to give 1,381 pixels globally representative of these five variables and their combinations. Forest cover and air temperature were not explicitly included in this subsetting, although a large range of both was covered. Figure 1 shows their location. 

All of the global products described above cover summer 2002 to summer 2008 and so this period was used for the analysis. Clifford (2010) showed that different SWE estimates had different spatial patterns, but the reasons for these were not obvious from the maps. In the present study an attempt was made to compare the products by binning mean monthly SWE by major river basins (Lena, Ob, Mackenzie and Yenisey). All that could be said from this analysis is that the three SWE products differ when aggregated up to the basin scale in terms of peak SWE (by as much as a factor of 2) and melt timing. The reasons for the differences and which is more accurate cannot be determined from these climatologies. These differences will be directly driven by the land surface and meteorological properties and so these are explored in the rest of this paper on a 1 o pixel basis rather than the spatial patterns, to more easily determine the mechanisms behind the differences.
section_name: Sites

section_content: Although a comprehensive evaluation of SCA is not the focus of this paper, an attempt was made to compare MOD10 and MODIS-SSM/I. If they performed equally well, the coarser resolution (and so smaller file) and cloud tolerant MODIS-SSM/I could be used to assess the SWE products; however, the blended product could suffer from the same issues as the SWE estimates described above. Time series of MOD10 and MODIS-SSM/I were plotted along with SSM/I SWE (as this is related to the microwave component of the blended product) and meteorological data, and snow start and end dates examined. Initially a visual comparison was made. 

There are many 500 m MOD10 pixels within a 1 o pixel, each of which can be independently cloud covered or masked for other reasons. Each 1 o pixel can thus have different numbers of MOD10 pixels contributing on different occasions; the greater the number the more representative the average SCA will be. The number of contributing MOD10 pixels, n, can be used as a measure of uncertainty, σ sca, using the following equation;

Where N is the total number of MOD10 pixels within the 1 o pixel.
section_name: Snow cover

section_content: MOD10 and MODIS-SSM/I snow start dates frequently disagreed. Of the 1,381 pixels and six years of overlap of the two products, 58% of snow start dates were successfully estimated (correct within sixteen days, two MODIS-SSM/I time steps). The failures were entirely due to missing MODIS-SSM/I data before snowfall. Of these, only 4. 8% agreed within eight days (a perfect agreement within the MODIS-SSM/I resolution) whereas 53% had MODIS-SSM/I predict a snow start date more than eight days before MOD10. Figure 2 shows two examples. In figure 3. 2. 1 there has been an early snowfall that completely melts (around the 12 th November) before a new snowfall that persists throughout winter. This melt and new fall occurs between MODIS-SSM/I measurements and so would give too early a snow start date. More worryingly, figure 3. 2. 1 shows a period of quite certain low SCA from MOD10 (late October to early November), whilst MODIS-SSM/I suggests snow cover should be complete. The Princeton air temperature was above 0 o C for both of these periods, making it likely that the snow would have melted, suggesting that MOD10 is correct. The snow end dates generally agreed better than the snow start date but did occasionally have MODIS-SSM/I drop to 0% before MOD10. Reasons for these differences are beyond the scope of this paper (perhaps ground frost or wet snow?). Therefore, if we trust MOD10 and Princeton air temperature (though the accuracy of Princeton temperature is not truly known (Sheffield et al. 2006 )), it would appear that MODIS-SSM/I will not always get the snow start dates right and cannot be relied upon for testing the SWE products. Only MOD10 will be used for the rest of this investigation, accepting that data will be missing due to clouds and the Arctic winter.
section_name: Snow cover results

section_content: Most global comparisons between SWE estimates have looked at global maps of monthly SWE and the snow dates (Clifford 2010, Blyth et al. 2010). Although these maps show the areas of disagreement, the reasons are not easy to determine. Time series of the three SWE variables, MOD10 SCA, air temperature and snowfall were first plotted to allow a visual investigation of their relative behaviour. There are three areas we should examine; the peak accumulation, vital for modelling the spring and summer hydrology as well as insulating the soil, the snow start and end dates, important for long term trends in snow's effect on surface albedo and growing season length, and the overall shape of the SWE seasonal shape, especially the melt period where the shape will determine whether water infiltrates the soil or runs off. This paper is most concerned with peak SWE and melt rates. Although snow dates are important for the surface albedo, they are the most difficult case for microwave observations, being shallower and occasionally wet.
section_name: SWE time series

section_content: Figure 3 compares the three SWE products to each other and meteorological data. Snow cover has been left off for clarity and will be addressed in section 3. 3. 4. As long as the air temperature remains below 273 K, SWE should be roughly equal to the cumulative snowfall. Some snow mass will be lost through melt and sublimation (at these scales redistribution of snow will be negligible) and these will depend on meteorology and land cover, but it gives a rough guide and an upper limit. These graphs are for the ABACUS site near Kevo in Finland (the 1 o centred on 69. 5 o N 27. 5 o E) and the BERMS (Boreal Ecosystem Research and Monitoring Sites) site in Canada (the 1 o centred on 58. 5 o N 104. 5 o ), both extensively studied areas for which we would expect the winter precipitation data to be accurate. They were chosen as illustrative of some common cases. This could well be due to the issues of the Chang algorithm in deep snow, as discussed in section 2. 1, although it occurs at a lower SWE than expected (Chang et al. 1987). 

At BERMS the SSM/I SWE is greater than all other estimates. It has been reported that all Chang based SWE estimates greater than 120 mm are spurious (Clifford 2010), possibly due to the forest correction or incorrect grain size. According to the BU-MODIS land cover map, Kevo is 0% forest covered whilst BERMS is 70%. At BERMS this will cause a doubling of the SSM/I estimate and the SSM/I estimate does seem to be double AMSR-E and Globsnow. 

Therefore it would appear that SSM/I and AMSR-E can suffer from saturation and that SSM/I is susceptible to exaggerated forest corrections whilst Globsnow appears to give similar peak accumulations to Princeton snowfall. Section 3. 5 will explore these effects in more detail.
section_name: SWE time-series results

section_content: The Chang based estimates often showed a significantly different shape to Globsnow, as shown in figure 3. The Chang based estimates show a sharp rise at winter's end that is not evident in Globsnow, whereas in figures 3. 3. 1 and 3. 3. 1 Princeton suggests there is little new snowfall at this time. Although in figure 3. 3. 1 Princeton does show new snowfall, this occurs when the daytime air temperature is above 273 K and thus when we would expect snow to be melting rather than accumulating. 

This was a common feature of the Chang based products, particularly for AMSR-E that showed it in the majority of areas examined. It seems to occur as air temperature rises above 273 K and then decreases again. As AMSR-E and SSM/I make measurements at night, the snow could well have refrozen. Such melt refreeze cycles are likely to dramatically alter the snow-pack structure, possibly creating ice lenses and thus increasing the effective grain size. This will be missed by SSM/I's static algorithm and could exceed AMSR-E's maximum grain size of 1 mm, thus leading to SWE overestimation by the Chang algorithm. 

We believe these spikes to be artefacts and Globsnow to provide the more accurate SWE seasonal shape, particularly during melt. These spikes substantially increase the peak annual SWE (though only during melt) and so would mask an underestimate of SWE and lead to a much more abrupt melt if ignored. The higher temporal resolution products show occasional spikes (particularly AMSR-E), but these can easily be removed by median filtering.
section_name: Seasonal shape

section_content: To determine the start and finish of the winter snow period (referred to as the durable snow cover by Kitaev et al. (2002) ) an algorithm was written to search for the first and last dates of permanent snow cover estimated by each snow product (SSM/I, AMSR-E, Globsnow and MOD10). First noise was removed using a five day weighted median filter (Brownrigg 1984 ). The weightings were set by the uncertainty. For MOD10, if the SCA uncertainty was less than 33%, weight three was allocated. For uncertainties between 33% and 66% weight two was allocated, otherwise weight one. For Globsnow (AMSR-E and SSM/I don't provide uncertainty estimates) if the SWE error was lass than 15 mm weight 3, between 15 mm and 35 mm weight two and greater than 35 mm weight one was allocated. 

The algorithm started from the 1 st February and tracked (forwards and backwards) until the SWE or SCA passed beneath a threshold (1 mm for SWE and 0. 5% for SCA), that point then being taken as the start or end date. If the signal didn't drop below this threshold within the year the retrieval was marked as unreliable and not used in subsequent analysis. This algorithm should avoid the small fall and melt events in autumn and spring but may be affected by noise. Only years with at least two weeks of snow cover (according to MOD10) around the 1 st February were used. 

Table 1 shows the mean error (bias) and root mean square error (RMSE) between the three SWE products and MOD10 (a negative bias is late, a positive early). All RMSE values were several weeks, showing that there was considerable spread in the differences, although we cannot be sure if this is due to inconsistencies in the SWE products or uncertainty in the MOD10 products along with weaknesses in the date-finding algorithm (and even in the definition of snow dates). SSM/I gave a much lower start bias than AMSR-E and Globsnow and a lower RMSE value, though a similar spread about the mean; however, the RMSE is still much higher than is desirable. We might expect AMSR-E and SSM/I roughly to agree, given their similar algorithms. That they don't suggests that the different temporal resolutions may be playing a role. SSM/I's eight day period will be less affected by noise than AMSR-E, Globsnow and MOD10's one day period. Often the Globsnow estimate did not start until after the first snowfall. These retrievals should have been flagged as unreliable but some will not have been (as the SWE occasionally started from 0 into the snow season); this may have increased the Globsnow errors. AMSR-E and SSM/I did not have these gaps. Takala et al. (2011) report that the snow start date is a particular weakness of the Globsnow product due to its reliance on the Chang algorithm to determine the start of the snow season. 

All products showed much better agreement for end date, both in bias and RMSE. The start and end dates in Globsnow will be determined by the Chang based snow start method of Hall et al. (2002) and the empirical snowmelt algorithm of Takala et al. (2009). Therefore we would expect the snow start date of Globsnow to be more conservative than SSM/I and AMSR-E, which it is. 

These results suggest that all three products struggle to determine the snow start date (though SSM/I gets it roughly right on average), but give a fairly good estimate of snow end date (though RMSEs are still two to three weeks). That all three products gave roughly the right snow end date however, does not necessarily mean they are giving an accurate measure of snow melt.
section_name: Snow dates

section_content: Whilst an inter-comparison of EO-derived snow products and meteorological re-analysis data will highlight any differences between datasets and can suggest reasons for these, we cannot know for certain which is correct without comparing against direct measurements. The three SWE products were compared to a point measurement of snow depth described in section 2. 7. The SWE products were used at their native 25 km resolution, although even then the difference in scales could cause disagreements that will be difficult to account for (Chang et al. 2005).
section_name: Ground data

section_content: Figure 4 shows the comparison of the three SWE products to ground measurements. Unfortunately SSM/I ends in early 2008 and so we cannot definitively evaluate its performance here. Figure 3, however, showed that in previous winters it behaved in a similar manner to AMSR-E over the ABACUS site near Kevo, and there is no reason to believe it would not do so for 2008 or 2009. 

This comparison suggests that AMSR-E is saturating around 80 mm and that the spike in latewinter/spring is spurious. The agreement between Globsnow and the ground data is good, especially given the assumption of constant density over time when converting from depth to SWE. All fluctuations shown by the ground data are mirrored by Globsnow, though not always of the same magnitude. The difference in fluctuation magnitude is not surprising given the difference in scales. 

Other than a spike in November 2008 (possibly due to scale differences, perhaps a passing snow drift) Globsnow gets the correct peak accumulation (within 10%) and melt date (within two days). This strongly supports the conclusions of section 3. 3, namely that AMSR-E and SSM/I saturate around 100 mm and show spurious spikes at the onset of melt, suggesting that Globsnow is the superior product. 

A miscalibration of the ground instrument is evident, giving negative snow depths in autumn; possibly due to seasonal senescence or bending of vegetation. This makes it difficult to more precisely evaluate the snow start date accuracy.
section_name: Ground data results

section_content: Once the general behaviour of the three products had been determined and potential issues highlighted, the mean SWEs for each month, year and product were calculated and compared (giving 72 separate SWE values per pixel, although summer values were uninteresting). The previous section suggests that, away from forests, SSM/I and AMSR-E estimates saturate around 100 mm SWE but may have a spurious spike at the onset of melt. Scatterplots were produced of each SWE product for each month to see if these features describe the products' global behaviour. Precipitation was not used in this analysis as it was considered unreliable away from intensively studied regions. Splitting the data into months avoided artefacts during the melt period (see section 3. 3. 3 and 3. 4. 1). Section 3. 3 suggested that the forest correction could be causing spurious SSM/I values and so the comparison was separated into for pixels with less than and greater than 10% forest cover (according to BU-MODIS). 

Figures 3 and 4 suggest that AMSR-E and SSM/I should saturate with respect to Globsnow over unforested areas; however, when scatterplots of the three products were produced for mid-winter (December-February) no such saturation was apparent (figure 5 ). At higher SWE values, AMSR-E and SSM/I diverged from Globsnow, but with no clear saturation of any product. Thus factors other than forest cover and saturation must be contributing to the differences. 

The comparison of AMSR-E to SSM/I (Figure 3. 5) seems to show a bimodal distribution rather than the more random spread when compared against Globsnow. Points with forest cover less than 10% are grouped around the 1 to 1 line, showing that AMSR-E and SSM/I agree. Points with cover above 10% have a shallower gradient, so SSM/I's estimate is higher than AMSR-E by a constant factor. This difference could well be due to the difference in forest corrections of the two products and the additional use of MOD44B in AMSR-E's derivation. Such a grouping of forested and non-forested points was not apparent when comparing AMSR-E and SSM/I to Globsnow (though forested pixels showed a very slightly higher SSM/I estimate). 

This analysis suggests that, whilst the forest correction does have an effect on SWE estimate differences, other factors are having an impact.
section_name: Globally representative analysis

section_content: Microwave estimates of SWE are sensitive to grain size (Chang et al. 1987). SSM/I assumes a constant value of 0. 35 mm, AMSR-E allows the value to vary between 0. 2 mm and 1 mm, based on surface temperature, Globsnow attempts to derive the value at each ground station and interpolate. These different approaches could cause large SWE differences. 

Snow grains are known to grow when there is a large thermal gradient across the pack (Slaughter and Crook 1973). Strong thermal gradients will cause large snow grains, of the order of one to several millimetres (Chang et al. 1982, Davenport et al. 2012), larger than is assumed in the derivation of AMSR-E and SSM/I. Clifford (2010), Brodzik et al. (2007) and Boone et al. (2006) cited these large grain snow-packs as a potential source of systematic bias in Chang based SWE estimates. To investigate this, differences between the monthly SWE estimates were compared to the thermal gradient. 

The temperature of the surface under the snowpack is needed to calculate the true thermal gradient. However, other than at a few scattered test sites, this is not available. Instead the air temperature above the snowpack was used (mean for a given year and month). This was compared to the difference between all three products each month for areas with 0% forest cover. This is equivalent to a continuous form of the Sturm snow classes (Sturm et al. 1995). 

Figure 6 shows the difference between SWE products against air temperature. The SSM/I -AMSR-E difference displayed no pattern and so is not shown. These two products are likely to be affected by the same processes and so this is not surprising. It can be seen that the Chang based methods roughly agree with Globsnow for air temperatures above -35 o C (though with some higher SSM/I estimates beneath -15 o C), but beneath these the Chang estimates tend to be larger than Globsnow. This is what we would expect if cold air temperatures were causing grain growth above the values used for SSM/I and AMSR-E's derivation. It is possible that the colder areas are more remote and so be further away from Globsnow ground stations, however, plotting the distance to the closest ground station against air temperature showed that this was not the case. There is considerable spread in the data, leading to a weak, statistically insignificant correlation between air temperature and difference in SWE estimates and this is likely to be due to the simplicity of the metric employed here (grain growth depends on the snow depth and air and soil temperature fluctuations throughout winter, which would require a full snow model to fully capture) and uncertainties in the land cover map used to mask forests (Herold et al. 2008). It does strongly suggest that very cold air temperatures (< -35 o C) cause the Chang based methods to overestimate relative to Globsnow. This could be due to processes other than grain growth, but Globsnow will account for them by altering the effective grain size (as a microwave scatterer). 

The site near Kevo had an average air temperature of -11 o C and BERMS -20 o C for December 2002-2008, both of which are high up the x-axis of figure 6 and so unlikely to show the overestimate of SWE by the Chang methods. It seems likely that figure 5 did not show saturation because of overestimates caused by grain size effects. 

The grain size maps derived as part of Globsnow are recorded and provided on request. These give the effective snow grain size (as used by the HUT model) for each 25 km pixel on each day. Here "effective" is used because the value could include other effects (liquid water, vegetation, soil temperature variations, etc. ) that act to change the relationship between SWE and the microwave signal, and thus may not have a strictly physical meaning. Whilst the effective grain size (as a microwave scatterer) may not be equivalent to the physical grain size, they might be expected to be related. To test this the mean grain size was calculated for each year and month and compared to the air temperature and difference between SWE products. 

Disappointingly the scatterplots of SWE difference against grain size according to Globsnow (not shown) did not show as strong a relationship as did those for the air temperature. The relationship between Globsnow grain size and air temperature was also not clear. It would appear that the effective grain size of Globsnow is not related to the same parameters that are the differences between SWE products and air temperature. The data used to derive this grain size are not readily available and so further investigation of this was not performed. For the rest of this paper Globsnow's grain size will only be taken as a parameter affecting Globsnow's derivation rather than a value with a strictly physical meaning. That is not to say that it doesn't have one, but that it is first and foremost a parameter fitted by the model and can be different from the actual grain size (Hall et al. 1986 ).
section_name: Snow grains

section_content: To determine the relationship between forest cover and SWE estimates the differences between the SWE values were compared to BU-MODIS forest cover. To avoid the potential influence of large grains, this was only performed for areas with warmer air temperatures. A threshold of -35 o C was chosen on the basis of the data shown in figure 6. Figure 7 shows the relationship between the difference between SSM/I and Globsnow against forest cover and the same for AMSR-E. There appears to be a weak relationship between the difference between SSM/I and Globsnow and forest cover up to 50% (the maximum correction factor value of 2). After this the difference drops again, suggesting that the forest cover is still having an effect upon the SWE estimate. There seems to be no clear relationship between the difference between AMSR-E and Globsnow and forest cover. Testing the goodness of fit of a straight and a sloped line to the differences between SSM/I (or AMSR-E) and Globsnow for forests with over 0% and less than or equal to 50% cover showed that it is not a statistically significant relationship. 

There is considerable scatter, making this relationship less significant than that for the thermal gradient, but it suggests either that Globsnow is not appropriately correcting for forests, or that the correction in equation 2 used for SSM/I is not appropriate, as proposed by Clifford (2010). Some of this scatter may come from uncertainties in the land cover maps used. Globsnow and the Chang based methods use different land surface maps and previous studies have shown that these can differ wildly, with either being wrong in certain locations (Herold et al. 2008).
section_name: Forest cover

section_content: The Chang algorithm, used by AMSR-E and SSM/I, was developed and tested over three sites, in Russia, Canada and the USA. The coordinates for these are given in Foster et al. (1980). These were chosen as relatively simple cases (Foster et al. 1980, Chang et al. 1982 ) and so may not be representative of the global snow-affected land areas. To test this a selection of surface and meteorological properties were extracted for the three areas and compared to the same values for the 1,381 1 o pixels used here. 

The properties tested were mean peak annual SWE, calculated from Globsnow; mean elevation and topographic variation (standard deviation of elevation within the site), calculated from ASTER G-DEM; forest cover, from BU-MODIS; and the mean January air temperature, from Princeton. Table 2 shows that the three sites used by Chang et al. (1987) are all much warmer than those used here (which should be representative of global snow-affected areas), have much shallower snow-packs and much lower forest covers. The topography alone seems not to be significantly different. It is thus not surprising that the same forest correction, saturation and large snow grain effects were not apparent during the algorithm's creation.
section_name: Chang algorithm calibration

section_content: The examination of Globsnow v1. 0 time series revealed occasional abrupt changes that could well be spurious. In an attempt to determine their origin these were compared to air temperature from Princeton (at that pixel and at contributing ground stations), effective snow grain size and which ground stations contributed. Princeton estimates of precipitation were not used as these could be unreliable in winter.
section_name: Globsnow artefacts

section_content: Figure 8 shows two locations with noticeable abrupt changes in Globsnow. The coloured bars along the x-axis indicate which ground stations (all within 500 km) were contributing each day. These same jumps are not apparent in AMSR-E or SSM/I and seem unlikely to be real. All the jumps appear to correspond to gaps in the nearest ground station's contribution. Interestingly, the effective grain size does not show a corresponding jump and occasionally shows a jump without much change in SWE, therefore it is likely that the SWE jumps are caused by a change in the interpolated ground snow depth map before it is combined with the EO data. The Globsnow values of the seven pixels within the 1 o pixel were examined and all showed the same jump, therefore it is not due to individual pixels dropping out. The temperature profiles of the pixel of interest and nearby ground stations are similar so it is unlikely to be due to different snow grain changes. Although it seems that most jumps are due to which ground stations contribute, they could also be due to changing weights given to the EO and ground estimates in the final combination, or else to temperature fluctuations altering the snowpack in a way that is not accounted for by the algorithm. Unfortunately the weightings are not recorded (Luojus pers comms) and so we cannot say for certain whether it is due to this or to temperature fluctuations; further investigation is needed. 

A full analysis was not attempted with versions 1. 2 or 1. 3 of Globsnow due to their recent release, but selected areas were examined to determine whether they too suffered from the same artefacts. A known issue has been introduced to v1. 2 over America due to a bug in the assimilation code and this should be solved in v1. 3 (Luojus pers comms) which has been released since this work was completed. 

Whilst v1. 2 did not show the same jumps as v1. 0 at the locations used in figure 8, figure 9 shows that it does suffer from similar artefacts. These do not necessarily occur in the same areas and times as in v1. 0, but they are similarly puzzling. A quick look at v1. 3 reveals that it also suffers from similar jumps. Therefore the newest version has not eliminated these artefacts and they should be looked at in more detail. Figure 9 shows that weekly v1. 2 product suffers from the same jumps as the daily product. The edges of the jumps were slightly smoothed, but the same dramatic, spurious shifts in SWE were evident. This smoothing may even mask the fact that these jumps are erroneous. These jumps were not common, only appearing in around 1% of hydrological years examined (1,381 pixels of Globsnow data, 1978 Globsnow data, -2008)), but they could cause issues when using Globsnow to validate model results, and especially if trying to assimilate Globsnow data into a model.
section_name: Globsnow artefacts results

section_content: This paper has compared three global SWE products, SSM/I, AMSR-E and Globsnow v1. 0, to each other, snow cover, ground measures of snow depth and meteorological data. The analysis was performed over 1,381 1 o pixels representative of the latitude, longitude, peak SWE, elevation and topographic variation of global snow-affected areas. 

The results suggest that the Chang algorithm based estimates (SSM/I and AMSR-E) can suffer from saturation above 80 mm -150 mm SWE, can show spurious spikes during melt (probably due to melt-refreeze), can show overestimates in the presence of cold air temperatures (likely to cause large snow grains) and in forested areas. In all cases the difference between the Chang based products and Globsnow behaved as we would expect the error to behave. This suggests that Globsnow is the superior product for peak accumulation and the seasonal pattern of SWE development. It would have been worrying if the newer, more advanced Globsnow product were not found to be superior. The sites used to develop and test the Chang algorithm were amongst the very warmest, had the shallowest snow and lowest forest cover of snow-affected areas; it is thus not surprising that the initial studies did not encounter these issues. 

The single point of ground data used suggested that the Chang based products were saturating and showing spurious spikes during melt, whereas Globsnow seemed to give a good estimate. This was over an unforested area with a relatively warm air temperature and so could not be used to test these effects. 

An attempt was made to assess snow start and end date accuracy. The results suggest that the start date for all three products has large errors (over one month), although SSM/I seems to agree with MOD10 on average. All products gave much more accurate snow end dates, with a mean bias of no more than a week, although the spread was still about a month. Therefore, it seems that none of the SWE products would be suitable for testing snow start dates of land surface models; this is not surprising given that passive microwave retrievals are insensitive to shallow or wet snow. 

Whilst Globsnow appears to be the superior product, it shows occasional spurious jumps. These seem to be due to gaps in the contribution from individual ground stations (an issue not affecting the pure EO products). These jumps only affected about 1% of the hydrological years examined, but they will need to be avoided or else taken into account if Globsnow is to be used to validate LSMs or be assimilated. No attempt was made here to investigate the stability of Globsnow at the crossover between sensors (1988, 2002 and 2011). 

Having determined that Globsnow is currently the best global SWE product we intend to use it test the snow processes within a land surface model, taking note of the issues above. This work suggests that the SSM/I and AMSR-E products are not suitable for such tests, especially during snowmelt and will even struggle to give the correct peak accumulation. Globsnow seems to have overcome these issues.
section_name: Conclusions

section_name: Acknowledgements

section_content: This work was funded by the NERC National Centre for Earth Observation (NCEO). 

Thanks to the Globsnow team, especially Kari Luojus, for providing useful comments and extra details of their product and to the five anonymous reviewers for their thorough and helpful comments.
section_name: Acknowledgements