The operational analysis of snow depth at the Canadian Meteorological Centre is described. The analysis makes use of forecasts of precipitation and analyses of screen-level temperature to estimate snowfall and snowmelt for a global domain, and assumes persistence of the mass of the snowpack between melting and/or snowfall events. In addition, wherever snow depth observations are available, these are incorporated using the method of statistical interpolation, performed every 6 h on a ⅓° grid. Correlations between two observations or between observations and grid points are taken as functions of spatial separation in both the horizontal and vertical with an e-folding distance of 120 km in the horizontal and 800 m in the vertical. Observations undergo three separate tests designed to eliminate 1) false reports of snow when none is present, 2) systematically understated snow depths, and 3) reports that violate temporal continuity. Verification of the analysis is presented using three sources of independent data. The analysis obtained when observations are withheld is compared with quality-controlled snow depth reports, with the result that the analysis shows more skill than climatology in all regions and periods examined. The analysis is also compared with Special Sensor Microwave/Imager snow cover, where it again shows more skill than climatology. Finally, a verification over a 134-week period using the NOAA weekly hemispheric snow cover product based on visible imagery is presented. Major differences between analysis and independent data are explained.
A need exists in the numerical weather prediction (NWP) community for real-time snow information, specifically, the geographic distribution, volume, and density of the snowpack. The current global snow-monitoring effort is discussed in Robinson et al. (1993). Real-time snow monitoring began in the late 1960s when the U.S. National Oceanic and Atmospheric Administration (NOAA) began producing weekly analyses of Northern Hemisphere continental snow cover from visible satellite observations. The digitized version of this analysis, hereafter referred to as the visible product, has been used in the surface energy budget computations of NWP models (Petersen and Hoke 1989). Model forecasts of surface air temperature, particularly the daily maximum temperature, are very sensitive to snow cover. However, the poor spatial and temporal resolution of the visible product has led to deficient forecasts in individual cases (Mitchell et al. 1993). More recently, progress has been made mapping snow cover using data from passive microwave radiometers flown aboard the Defense Meteorological Satellite Program (DMSP) satellites (Grody 1991). This product represents a substantial improvement in spatial and temporal resolution over the digitized weekly snow cover analysis of NOAA.
Useful though it may be for many applications, snow cover information is inadequate for the computation of surface radiation fluxes in NWP models. In this context, it is desirable to know the depth and density of snow at the model grid points. Snow depth is needed to refine estimates of the surface albedo. Robinson and Kukla (1985) report that albedo increases rapidly with snow depth, up to about 15 cm. At the Canadian Meteorological Centre (CMC), snow depth is combined with vegetation information to estimate albedo for a global domain. Mailhot et al. (1997) use snow depth to partition available heat between melting snow and sensible heat fluxes. If the density of the snow is available along with the depth, one may more accurately parameterize the melting process in hydrological computations. In addition to the greater utility of snow depth, the choice of snow depth rather than snow cover as the analysis variable offers a compelling advantage. With snow depth and density information, it is possible to exploit temporal continuity by applying conservation of mass from one analysis time to the next. This is not possible if only the snow cover is known.
The literature on the subject of real-time snow depth analyses is somewhat sparse. The U.S. Air Force daily global snow depth analysis has been documented by Hall (1986) and later refined by Armstrong and Hardman (1991). It makes use of the available observations of snow depth and passive microwave satellite data; however, where neither of these data sources are available, the analysis trends toward climatology. Since snow depth anomalies from climatology are frequently quite large, it is desirable to find an alternative to climatology to fill data voids. One possibility is to use information from an atmospheric model. Buchhold (1996) describes how this is done at Deutscher Wetterdienst (Offenbach). There, snow depth is a prognostic variable of the global NWP model, and the operational global snow depth analysis incorporates the forecast snow depth in data-sparse areas.
This paper describes the global objective analysis of snow depth that is produced operationally at the CMC. The analysis has been designed to address the issues discussed above: it has a spatial and temporal resolution comparable to that of the passive microwave observations, it provides the required information on snow depth and density as well as snow cover, and it does not rely on climatological information. The next section provides details of the method.
Description of the snow depth analysis
All snow depths received from the Global Telecommunication System in SYNOP (land surface synoptic report) format are used in the analysis. It is the practice of some North American stations to indicate snow depth only in the remarks section of a METAR (aviation routine weather report), according to a prescribed format. Where available, snow depths are decoded from the METAR reports to supplement the SYNOP dataset. Observations from a 6-h period centered on the analysis time are first screened to select those with a snow depth measurement. If a station has reported snow depth more than once over the 6-h period, only the report closest to the synoptic hour is used in the analysis. The World Meteorological Organization (WMO) code allows for snow cover that is patchy in nature to be reported. The analysis decodes this as 1 cm, the smallest nonzero snow depth allowed by the code.
A map of stations that reported a nonzero snow depth from 1 October 1996 to 31 May 1997 and whose report was accepted by the quality control is shown in Fig. 1. Many stations report snow depth only once per day. This fact explains the daily variation in the number of reports from a minimum of approximately 200 stations at 1800 UTC to a maximum of about 600 stations at 0600 UTC (during boreal winter). Stations cease to report snow depth if no snow is present, resulting in a minimum in the number of reports in the summer. For the Southern Hemisphere, this seasonal paucity of data is clearly seen in Fig. 1, although it should be pointed out that stations reporting snow depth are few even during austral winter.
The analysis is computed on a global, ⅓° latitude–longitude grid. Care must be taken to ensure that none of the observations used are inconsistent with the mean topography at this resolution. Figure 2 is a histogram of station elevation minus mean topography height for all the stations of Fig. 1. The topographical data are the National Geophysical Data Center (1988) 5-min elevations, which were area-averaged onto the analysis grid. The reference station heights required were obtained from the index of station data (WMO 1996). It is noteworthy that a substantial negative bias exists in this distribution; that is, many more stations are below the mean height of the topography than are above it. This raises the question of representativeness of the observations. The analyzed snow depth should represent conditions at the mean topography elevation for the grid square. This will not be the case, however, if the analysis uses observations from stations whose elevations are significantly different from that of the mean topography. To resolve this problem, the data from stations more than 400 m above or below the mean topography height are not used. Even after eliminating these large differences, however, the median difference of the remaining stations is −26 m, indicating that there is a low-elevation bias in the station data.
Prior to each analysis, a background snow depth is computed. The most recent (i.e., six hours old) snow depth field is the basis of the background field computation. This is an appropriate choice since it can be shown from quality-controlled snow depth observations that the average 24-h lag correlation for stations during Northern Hemisphere winter (December–February) exceeds 0.85. Thus, snow depth is a conservative quantity, and this property is exploited here.
In a complete model of the snowpack, one should take into account precipitation of all types, sublimation/condensation, snowmelt, and the metamorphosis of the snowpack (see, e.g., Lynch-Stieglitz 1994). Here, an operational forecast of precipitation accumulated over the 6-h period since the previous analysis time is used to estimate snowfall. The results presented in sections 3 and 4 were produced when the operational global NWP model described in Ritchie and Beaudoin (1994) was used. It had a resolution of 0.9°, and the data assimilation scheme used to specify initial conditions for the model was as described in Mitchell et al. (1996). The precipitation is added to the most recent snow depth field if the analyzed screen-level temperature is zero or less. The use of precipitation at a resolution of 0.9° means that where the snow depth varies on scales smaller than this, the analysis may be deficient. This is particularly true for the “lake-effect” snowfalls occurring downwind of the Great Lakes, where the station density is inadequate to correct for the low-resolution precipitation information employed.
In addition to the snow depth, which represents an average over each grid square, the analysis scheme estimates the mean density of the nival layer at each grid point. Together, these two parameters may be used to compute the snow water equivalent as required. Fresh snow settles and undergoes metamorphism due to a variety of processes (Sommerfeld and LaChapelle 1970), resulting in a gradual increase in density over time. Accordingly, the density ρ is updated every 6 h using the formulation of Verseghy (1991)
where ρ−6 is the most recent density, ρa is the asymptotic value of the density due to aging, and the e-folding time τ is 100 h. Fresh snowfall is assumed to have a density of 100 kg m−3. The ρa of 300 kg m−3 used by Verseghy is appropriate for most terrain and vegetation types. However, the empirical data of McKay and Gray (1981) clearly show that the mean density of the snowpack in the boreal forest remains less than about 210 kg m−3 until the start of the spring thaw. Accordingly, the vegetation index of Wilson and Henderson-Sellers (1985) is consulted so that at grid points where the vegetation is primarily of the needleleaf variety, ρa is set to 210 kg m−3 and elsewhere it is set to 300 kg m−3.
Following every snowfall, the updated density is calculated using a mass weighting of the density of the previous snowpack and that of the new snow. When the temperature is above 0°C, the density is increased at the rate of 0.5 kg m−3 K−1 h−1 to account for the increase in density associated with meltwater infiltrating the snowpack and refreezing. If the density exceeds ρa, the snow-aging algorithm (1) is not applied, but the density may still increase (decrease) if there are further snowmelt (snowfall) events. However, the density is limited to a maximum of 550 kg m−3.
A simple melting algorithm is applied when the analyzed temperature is greater than zero. Mass is removed from the snowpack at a rate that varies linearly with temperature. The constant mass removal rate, expressed in terms of water equivalent, is 0.15 mm h−1 K−1. Partly due to a lack of appropriate real-time information, the current snowpack model has been developed using a number of simplifying assumptions. In particular, it should be noted that the precipitation algorithm ignores liquid and freezing precipitation as well as condensation; the melting algorithm neglects solar radiation, heat fluxes from the ground, and losses due to sublimation.
The analysis is updated every 6 h using the method of statistical interpolation (see Daley 1991, chapter 4; Lorenc 1981). In the first step, the background field described in the previous section is interpolated horizontally to the observation locations using a bilinear interpolation scheme, and the observed increments from the background, fi, are calculated at every observation location i. Using a carat to denote an analyzed quantity, the analyzed increments, f̂g, at each grid point g are then computed using the linear-weighted sum
where wi are the optimum weights. The vector of weights w is given by
where P is the correlation coefficient matrix of background field errors between all pairs of observations, O is the covariance matrix of observational errors (normalized by the background error variance) between all pairs of observations, and q is the vector of correlation coefficients μig of the background field error between the observations and the grid point.
The correlation coefficients of P and q are vital to the method and are assumed to have the form
where rij is the horizontal separation of points i and j, and Δzij is the vertical separation of the points. The horizontal correlation function is the second-order autoregressive function (Daley 1991, chapter 4):
with c = 0.018 km−1, corresponding to an e-folding distance of 120 km. The function is plotted in Fig. 3. The form of the vertical correlation is
with the vertical length scale h set to 800 m.
Figure 4 illustrates the effect of the vertical correlation on the analyzed response to two mountain stations. The horizontal correlation, since it is isotropic, consists of concentric circles centered on the stations (not shown). The topography, shown in Fig. 4a, modulates the response, restricting the influence of the station in the Sierra Nevada of California to higher elevations to the east, and also preventing the station in western Montana from influencing the Rocky Mountains to the east of the station.
To solve (3) one must specify the observation error covariance matrix. It is assumed that errors in neighboring observations are uncorrelated, in which case O reduces to the identity matrix multiplied by o2/b2, the ratio of the observation error variance to the background error variance. Observation error includes observer error, instrument error, and the errors of representativeness due to the variability of the snow depth field on scales that cannot be resolved by the present analysis scheme. For a given station, the latter two error types can be expected to be temporally correlated, introducing local biases in the analyzed snow depth field. Since the snow depth analysis is the basis of the background for the subsequent analysis, in general, the background error cannot be assumed to be independent of the observation error. This fact makes it impossible to estimate the variance ratio using the traditional approach (see, e.g., Mitchell et al. 1990). Accordingly, the variances were chosen arbitrarily to be 10 cm2 for the background and 6 cm2 for the observations. The specification of variances, which determine the relative weight given to the observations with respect to the background, is not as important as in atmospheric data assimilation where the background is usually assumed to be independent of the observations. Here, the use of persistence in the computation of the background field means that, over time, any station reporting regularly will come to dominate the analysis locally regardless of the values of the variances.
Snow depth is an inhomogeneous field that varies with elevation, slope and aspect of the terrain, vegetative cover, exposure to wind and direct sunlight, and a host of meteorological factors, such as air and ground temperature, moisture, and precipitation on all scales. This lack of homogeneity makes the quality assurance of observations extremely difficult. Variability on small scales limits the usefulness of neighboring reports when checking for errors. In some circumstances, it is possible to use ancillary information to identify gross errors in the data, but in general one must be content to reject only those reports that flagrantly violate spatial and temporal continuity.
With the advent of completely automated meteorological observing stations, there has been a discernible increase in the occurrence of gross errors in the real-time reports of snow depth. A preliminary check is made to attempt to detect the most prevalent of these errors, the reporting of snow where none is present. If one assumes that stations with snow cover report snow depth at least once per day, then an appropriate consistency check would be to verify the reported 24-h precipitation (if available) and screen-level temperature over the 24-h period. Hence, if the a priori information (i.e., background) indicates there is no snow, then the reported snow depth is rejected if either the reported 24-h precipitation is zero or the minimum of the 6-hourly temperatures over the period is greater than 2°C. The 6-hourly temperatures are those of an objective analysis of screen-level temperature, which includes quality control of the temperature observations. The temperature analysis is described in section 2e.
Another useful source of information for checking observations is the snow depth estimate described in section 2b, which is essentially based on forecast precipitation and analyzed temperature. The error characteristics of this estimate are presented in detail in section 4. On average, the estimate systematically underestimates the snow depth. In addition, the error distribution (observed minus estimated) is skewed with more errors to the right than to the left of the mode. It is rare for an observation to be significantly less than the estimate. This fact provides an opportunity to remove outliers in this category. Therefore, observations that are less than the estimate by 40 cm or more are identified as erroneous.
There is also a final check applied to each datum using the statistical interpolation methodology. Following Lorenc (1981), an observation is rejected if it satisfies
where f̂i is the analyzed snow depth valid at the location of this observation and obtained by using only neighboring observations (if any), T is a user-specified tolerance (currently set to 5), and a2 is the theoretical analysis error variance given by a2 = b2(1 − Σ μiwi). Due to the small decorrelation length scales, this test often amounts to a check against the background, and since the background reflects previous reports from the station as described above, a check against the background is essentially a check of the temporal continuity of reports from the station. For an isolated station [a2 = b2 in (7)], the reported snow depth must be within 20 cm of the background to be accepted. Together, the three quality control checks described here reject an average of 3.7% of reported snow depths.
Screen-level temperature analysis
The observations used in the screen-level temperature analysis come from all available SYNOPs, supplemented by METARs if no SYNOP is available. The analysis is produced every 6 h on a global 0.9° grid. Observations from stations more than 400 m from the mean topography height (using topography averaged onto the 0.9° grid) are discarded to obtain a reasonably representative dataset.
The statistical interpolation methodology described in section 2c is employed for screen-level temperature as well. The background is the 6-h forecast of 2-m temperature from the global NWP model (Ritchie and Beaudoin 1994). The horizontal correlation function for temperature is a sum of third-order autoregressive functions given by
The vertical correlation is the same as that used for snow depth given in (6). The three important parameters of the method, b2, o2, and c are based on the 1000-hPa temperature statistics derived by Mitchell et al. (1990). The width parameter c and the background error variance b2 vary with latitude. The dashed line in Fig. 3 shows the form of the function (8) at 50°N, where the value of c is 0.008 km−1, corresponding to an e-folding distance of 425 km. The observations are checked using the test (7) above.
The resulting analysis represents the temperature at the mean level of the topography on a 0.9° grid. Prior to employing this analysis in the snow depth scheme, the analyzed temperatures are adjusted to account for differences in elevation between the 0.9° topography and the ⅓° topography of the snow depth grid. The adjustment assumes a lapse rate of 0.006 K m−1, which is the saturated adiabatic lapse rate for a temperature of 0°C and a pressure of about 800 hPa. Finally, the resulting 6-hourly temperature fields are linearly interpolated in time to yield hourly temperatures for use in the snowfall and snowmelt algorithms.
Spatial and temporal variation of the analysis
Figure 5 shows the analyzed snow depth for 0000 UTC 1 January 1997. The 2-, 30-, and 100-cm contours are shown with progressively darker shading for larger snow depths. Several relative maxima are indicated. These tend to occur at high elevations along storm tracks. The spatial variability is greatest where there are abrupt changes in elevation such as over western North America or central Asia. The surprisingly large snow depths over the Coast Mountains of British Columbia and the Rocky Mountains of Idaho resulted from colder than normal conditions over the latter two months of 1996 combined with a very significant snowfall event in late December.
Much may be learned about the characteristics of the analysis by comparing observations with an estimate based solely on precipitation forecasts and analyzed temperatures (see section 2b). This field, hereafter called the estimate, is equivalent to an analysis produced without benefit of snow depth observations. In Fig. 6, reported snow depths from two stations (dots) are compared with the analysis (solid curve) and the estimate (dashed line) for the period October 1996–April 1997. Station 71896 (Fig. 6a) is a Canadian station located just west of the Rocky Mountains, and 71943 (Fig. 6b) is located just east of the Rocky Mountains. The station locations are indicated by dots in Fig. 4a. Both stations are located in the boreal forest region, and so the mean density of the snowpack, shown in the upper panels, tends to be near the asymptotic value of 210 kg m−3 through the early winter. For 71896, the estimate agrees quite well with the observations until about the beginning of February, when a negative bias sets in and gradually worsens until the end of the season. By contrast, 71943 reports much more snow than the estimate right from the start of the season, and the bias in the estimate increases with time during the winter. The two stations illustrate two regimes of forecast error for precipitation. In the westerlies, the NWP model’s precipitation forecasts have a small systematic error for the western (windward) slopes of mountain ranges, but the systematic error is much more significant to the east (lee) of mountain ranges.
Figure 7 illustrates the operation of the quality control algorithm that eliminates reports that significantly understate the snow depth, as discussed in section 2d. The figure shows the reported snow depth from a station in eastern Canada (71728, at 48.3°N, 72.2°W) as well as the analysis (solid line) and the estimate (dashed line). The dots in the figure indicate that the snowpack was never reported to be more than 32 cm deep. These observations are difficult to reconcile with those from a station just 90 km to the east, which reported a maximum snow depth of 134 cm on 23 February 1997. All available information indicates the data are too low, with the apparent error increasing as the season progresses, then decreasing in the spring. The analysis follows the observations closely until 15 January. All subsequent reports are rejected by the quality control algorithm, freeing the analysis to follow the trend of the estimate with some influence from the nearby station. An investigation of this case revealed that snow depth reports from the station have been biased since 1992 when the snow course was moved from a location near a wooded area to a more exposed area adjacent to an airport runway, suggesting wind redistribution as the cause of the problem.
Extreme events, when they occur, test the responsiveness of an objective analysis. The heavy snowfall over Massachusetts on 1 April 1997 was such an event. Figure 8 portrays the time evolution of the analyzed snow depth for 28 grid points along 42.5°N. At about 72°W, the analysis goes from an absence of snow to a snow depth of 33 cm in just 18 h, in response to a METAR from Chicopee Falls (42.2°N, 72.5°W, reporting 34 cm). This is a rapid change, but the 6-hourly updates allow the analysis to respond appropriately with a smooth transition to the maximum value and a smooth return to zero during the melting phase. In this case, the observations are essential to the quality of the analysis since the maximum of the estimate was only 13 cm. The latter result was likely due to the fact that much of the snow fell when temperatures were above 0°C.
Comparison of the estimate with independent snow depth observations
It is of interest to evaluate the error of the snow depth analysis using independent data. Unfortunately, there is little reliable independent information on snow depth. Therefore, in this section, the error of the estimate is presented with the independent data provided by the quality-controlled snow depth reports actually used in the analysis. The error of the estimate is useful to gauge analysis error in data-void regions. Since almost all stations reporting snow depth are in the Northern Hemisphere, it is appropriate to examine the error characteristics during boreal winter.
Table 1 shows the errors of the estimate for the periods 1 October 1995–31 May 1996 and 1 October 1996–31 May 1997. The statistics for the two periods are averaged, and Table 1 presents the rms errors and biases by month for the Northern Hemisphere as well as 16-month averages for several regions. The regions are all bounded by the equator and 90°N, with the following definitions for the longitudinal boundaries: Europe, 0°–60°E; Asia, 60°E–180°; western North America (WNA), 100°W–180°;and eastern North America (ENA), 50°–100°W. As a benchmark, the U.S. Air Force global snow depth climatology (Foster and Davey 1988) was also verified against the same observations and the results are included in Table 1. For this comparison, the 12 monthly climatologies were interpolated linearly in time and bilinearly in space to yield a climatological snow depth referenced to the time and location of each observation.
The rms errors of Table 1 show the estimate to be more accurate than climatology for every month and region examined, although the difference is not significant at the 95% level of confidence for western North America, where the estimate has a large negative bias. In fact, averaged over the Northern Hemisphere, the estimate has a 22% smaller error than climatology. It is reassuring that the basic requirement for an objective analysis (i.e., that it show more skill than climatology) is satisfied even without assimilating any observations. On the other hand, the biases show that the estimate suffers from a systematic underestimation of the snow depth with the overall mean of −3.0 cm for the Northern Hemisphere. It is likely that most of this bias can be attributed to the underestimation of precipitation by the NWP model, which appears to be most severe over the mountainous western North American region and over Europe. However, another factor to consider is the assumption that the snow density is constant with depth and the associated mass-weighting of the densities when there is new snow. The effect of this assumption is that the larger the depth of the snowpack, the smaller the increase in depth for a given mass of new snow. In reality, the change in total depth due to a given quantity of new snow should be independent of the existing snowpack, but accomplishing this numerically requires treating layers of the snowpack independently and greatly complicating the real-time analysis scheme. Biases introduced by this limitation are corrected by the analysis where observations are available; elsewhere, the bias diminishes as the fresh snow ages.
Comparison of the analysis with SSM/I snow cover
The mapping of snow cover using microwave radiances measured by the Special Sensor Microwave/Imager (SSM/I) instrument has recently attracted attention as an information source to supplement the NOAA visible product. The original SSM/I snow cover algorithm of Grody (1991) has recently been refined by Grody and Basist (1996), and there has been a substantial effort to document the differences between the SSM/I analyses and the visible product (Basist et al. 1996). In this section, the Grody and Basist (1996) algorithm is applied to derive independent snow cover data against which the analysis described here is verified.
The verification covers the Northern Hemisphere for the period 1 November 1996 to 30 April 1997. The SSM/I-derived snow cover data (hereafter referred to simply as SSM/I snow) are referenced to a 384 × 384 polar stereographic grid (approximately 48-km resolution). There were two operational DMSP satellites during the period, F10 and F13, with equatorial crossing times of approximately 0900 and 0600 LST, respectively. Grody and Basist (1996) caution that to minimize a problem with detection of melting snow, only the morning satellite observations should be used. Since the morning passes of F10 are typically after sunrise, only the morning passes from the F13 satellite are used here. In order to compare data valid at the same time, no attempt was made to use SSM/I snow from previous days to fill orbital gaps. The elimination of all F10 data and evening passes from F13 resulted in a data frequency of approximately one observation every two days at 50°N, with lesser frequency of coverage farther south and greater frequency farther north. The verification was not carried out at grid squares where land accounted for 50% or less of the area. To compare with the SSM/I snow, the snow depth analysis from 1200 UTC each day was interpolated to the reference polar stereographic grid using a bilinear interpolation scheme, and a threshold of 1 cm was used to identify snow-covered grid squares. As in Table 1, climatology was verified against the same data. In the case of the climatology, it was not necessary to interpolate spatially since the grid of the U.S. Air Force global snow depth climatology is the same as the verification grid.
Table 2 presents the results of the comparison of the analysis and climatology with the SSM/I snow cover. Defining n2 as the number of cases where the analysis (climatology) showed snow cover and SSM/I did not, n3 as the number of cases where SSM/I showed snow cover and the analysis (climatology) did not, and n4 as the number of cases where both SSM/I and the analysis (climatology) showed snow cover, then the probability of detection (POD) is n4/(n3+n4) and the false alarm ratio (FAR) is n2/(n2+n4).
Table 2 shows that the analysis has a higher correlation with the independent data than climatology for every region and period examined. The analysis POD for the Northern Hemisphere is 0.94, meaning only 6% of snow cover was undetected by the analysis. Also, the analysis POD exceeds 0.9 for all regions except Europe. The many mountain ranges of Europe contributed to the low POD for the region. As defined above, Europe includes the Caucasus Mountains, the Zagros Mountains, the Taurus of southern Turkey, and the mountains of the former Yugoslavia. These four mountainous regions were associated with low PODs. It is curious that the region with the best POD, eastern North America, was also the region with the poorest correlation. In addition, the region with the best correlation, western North America, was also the region that had the poorest snow depth verification results in Table 1.
To help interpret the statistics for North America, Figs. 9, 10, and 11 show the snow cover frequency for the analysis and SSM/I for three periods during the winter of 1996–97. To facilitate the comparison between analysis and SSM/I, frequencies were based only on cases for which the SSM/I snow was available. For November and December (Fig. 9), the most apparent discrepancy is over eastern North America, where the analysis shows snow cover for more than 75% of cases for most of Quebec, Ontario, and Manitoba, whereas the SSM/I snow cover frequency is significantly less than this. Grody and Basist (1996) also noted this problem in a discussion of a mid-November case, attributing it to very shallow snow (less than 1 in.), which is below the sensitivity of the SSM/I instrument. Another possibility, cited by Robinson (1993), is that unfrozen soil beneath the snowpack contributes to the SSM/I underestimate during autumn. This type of error results in an elevated FAR for November, as seen in Table 2. For western North America, the two approaches agree quite well for regions north of 55°N, as well as for the Rocky Mountains, the Wasatch Range, the Sierra Nevada, and the Cascade Range. This relatively good agreement for snow cover in mountainous regions persists through the winter to late April and leads to higher correlations for western North America in Table 2. This result may appear to contradict the large rms error and bias for this region from Table 1, but it should be remembered that the in situ data used in Table 1 are not uniformly distributed through the region, and thus error estimates based on these data may not be representative.
The most striking feature of the January–February 1997 comparison (Fig. 10) is the difference in the region from the Great Lakes eastward to the Atlantic coast. In this region, the analysis shows snow cover frequencies of greater than 75%, while SSM/I shows significantly less than this. In fact, the region had snow cover throughout the period. Various researchers have reported that snow cover for dense coniferous forests is underestimated by passive microwave sensors (Hall et al. 1982; Basist et al. 1996; Chang et al. 1996), which may explain the difference over part of the region. Also of interest in Fig. 10 is the questionable lack of snow cover shown by SSM/I for coastal British Columbia, Alaska, Newfoundland, Labrador, and along the east coast of Greenland.
The comparison for March–April 1997 (Fig. 11) shows good agreement of the 0% contour running from western Oklahoma to Lake Michigan. The underestimation of snow cover east of the Great Lakes continues to be a problem for the SSM/I algorithm, as is the underestimation near some coastlines. Along the north shore of Lake Superior, for example, frequencies were as low as 4%, while just one grid length away, frequencies of 70% or higher were observed. Low frequencies near coastlines could be due to mixed-pixel (water–land) effects.
Comparison of the analysis with snow cover derived from visible imagery
The weekly NOAA visible product is a second independent source of snow cover information. In this section, analyzed snow cover will be compared directly to the visible product. The comparison spans 134 weeks from 2 October 1995 to 26 April 1998.
Two issues complicate the comparison. The first is the asynoptic nature of the visible product, which is based on satellite imagery from a 7-day period with preference given to the most recent image with clear-sky conditions. This means that the values at different grid points do not necessarily represent conditions on the same day. To make a comparison, therefore, one must choose the most appropriate analysis from the 7-day period. It was found by experimentation that an analysis from the seventh day of the period agreed best with the visible product.
The second issue is the question of how to transfer the snow depth analysis from the ⅓° latitude–longitude grid to the 89 × 89 grid of the visible product (approximately 190-km resolution). Since the visible product is hand drawn and then manually mapped to the 89 × 89 grid, it seems reasonable to adopt the same criterion used by the analyst. Essentially, a grid square is snow covered if 50% or more of its area is snow covered on the original map. When estimating areas, the analyst must also consider the following guidelines from the standard operating procedures (T. Baldwin 1998, personal communication): “The purpose is to show areas of snow cover: (the analyst should) include a box as snow if the immediate area totals one box, even if no individual box equals 50% coverage.” Therefore, mapping the analysis to the 89 × 89 grid using a threshold of 50% does not emulate the approach taken to digitize the visible product in every case. To allow for this difference, a grid square is defined as snow covered if 30% or more of its area has 1 cm or more of snow.
Figure 12a shows the correlations between the analysis and the visible product at each grid point of the 89 × 89 grid. Where there is no variability of the snow cover, the correlation is zero. The maximum correlation occurs at seven points where only one or two snow cover events were observed over the 134-week period and the products were in agreement. The correlations over most of northern Eurasia and much of northern North America exceed 0.8 and often exceed 0.9, indicating good agreement.
Figure 12b shows probabilities of detection of snow cover. Areas where snow cover was not observed during the period are hatched. Values exceed 0.9 for the majority of grid points where snow was observed. Over western Europe and the western United States, however, the probabilities are lower, resulting in lower correlations for these regions as well. Over western Europe, the transient nature of snow cover appears to make its detection by the analysis more difficult and could also lead to inconsistencies with the visible product if the latter is a day or more out of date. Over the western United States, it appears that snow cover at higher elevations is sometimes undetected by the analysis. There is evidence of this elsewhere in the figure as well (see, e.g., the European Alps, the highlands of Turkey, and the southern slopes of the Himalayas). The most likely explanation for this is that the analysis underestimates the extent of snow cover over highlands. The cause of the underestimation appears to be the temperature field (see section 2e), which corresponds closely to the NWP model forecast over data-sparse or data-void mountain regions. These model forecasts are too warm at high elevations, probably due to an error in the specification of the soil temperature. Work has begun to correct the problem.
One area where Fig. 12 shows a more serious problem, however, is over the Tibetan Plateau, where correlations vary between 0 and 0.3. Interestingly, the SSM/I snow cover algorithms also encounter difficulty in the region (Robinson and Spies 1994; Basist et al. 1996), although the cause is quite different. In the case of the analysis described here, the low correlations result from an excess of precipitation forecast by the NWP model over Tibet. Winter is the dry season for the region, with most of the plateau receiving less than 2 mm of precipitation on average during January (Domros and Gongbing 1988) as compared to model forecasts in the range 20–30 mm. The overestimation of precipitation over the Tibetan Plateau is a problem from November to March. The causes of this error are under investigation, but its effect on the analysis is to create a winter snowpack that persists until summer. The NOAA visible product, in contrast, shows that snow-free ground is typical in the region during the arid winter months. Unfortunately, no in situ observations of snow depth are available from the region to correct the analysis error.
Table 3 shows the probabilities of detection, false alarm ratios, and correlations for the snow depth analysis and climatology compared to the visible product. In Table 3, Northern Hemisphere refers only to land areas within the domain of the 89 × 89 grid (see Fig. 12). Also, grid points over Tibet have been excluded in order to give statistics representative of the remainder of the hemisphere. The results in Table 3 show that the analysis is better correlated with the visible product than climatology for all regions examined. As was seen above when SSM/I data were used, the analysis POD is lowest over Europe and highest in eastern North America. Snow cover undetected by the analysis over the grid domain (excluding Tibet) is 9% of the total, which is slightly higher than the 6% from Table 2. It is unclear how much of the undetected snow was due to the change of grid resolution since analysis errors were frequently found to be in areas where the analysis had some snow cover, but the total area was less than the threshold of 30%.
Figure 13 shows the area of snow-covered land over the grid domain expressed as a percentage of the total area with the Tibetan Plateau excluded. There is a tendency for the analysis to show more snow cover than the visible product during fall and less in spring. During fall, the analysis leads the visible product by 3 or 4 days. This difference could be the result of more frequent cloudiness necessitating the use of imagery that is, on average, a few days older than is commonly used during the rest of the year. In May and June, however, the analysis leads the visible product by about 2 weeks. This could be caused by the accumulated deficit in winter precipitation (see section 4a), resulting in the snow cover disappearing early. Also, the problem with the temperature field over highlands discussed above is likely to contribute to the error in spring since temperatures that are too warm will cause the analyzed snow to disappear too soon. During summer and winter, the differences in total snow-covered area are typically less than 3% of the total land area.
Figure 14 shows the weekly correlation between analysis and visible product (open circles) and between climatology and the visible product (closed circles). The correlations are computed over land points of the 89 × 89 grid with Tibet excluded. The analysis correlations are bounded by 0.60 and 0.92. The analysis shows more skill than climatology for 116 out of 134 weeks and performs best compared to climatology during summer and fall. The fact that the analysis typically has correlations near 0.9 during winter is encouraging.
Discussion and conclusions
A global, real-time analysis of snow depth has been presented. The method uses two sources of information: observations of snow depth from the synoptic observing network and snowfall estimated from NWP model forecasts of precipitation. The information in the observed snow depths is incorporated using statistical interpolation methodology. A background field for the statistical interpolation step is essentially obtained by combining a precipitation forecast with an analysis of screen-level temperature to estimate snowfall and adding the latter to the previous snow depth analysis. A temperature-based melting algorithm is also included. When there is neither snowfall nor melting, it is assumed that the mass of the snowpack is conserved.
The method is able to produce reasonable snow depths even for regions lacking in observations. Section 4a presents a verification of the snow depth analysis with observations withheld. It shows that rms errors of the estimate are better than those of climatology for every region and every month examined, and the error is 22% smaller than that of climatology for the Northern Hemisphere as a whole for the 1995/96 and 1996/97 snow seasons. However, it is clear from the verification results that the global NWP model suffers from a precipitation deficit that varies in severity with region. In addition, the bias estimated with reported snow depths may well be understated since stations are usually located in urban areas and/or in valleys (see section 2a), both of which diminish the representativeness of the data. The problem with biased and inaccurate precipitation forecasts is one of the great challenges of NWP research at present. Analysis quality can be expected to benefit with every future improvement in forecasts of winter precipitation.
The analysis scheme described here produces both snow depth and snow density fields. Since there are no real-time observations of snow density, an estimate based primarily on the age of the snowpack is made. The objective analysis step (section 2c) corrects the snow depth using the data but leaves the density estimate unchanged. Thus, it is implicitly assumed that the error in the background is due entirely to error in the snowfall accumulation algorithm, primarily the forecast of precipitation. Clearly, errors will be due as well to deficiencies in the melt and density algorithms and a host of physical effects not considered here. But it should be pointed out that the underestimation of precipitation by the NWP model also contributes to density estimates that are too large throughout the snow season. Without doubt, the single most important component of background field error is the bias inherent in the precipitation forecasts. It is tempting to try to reduce this systematic error by combining reported precipitation amounts with NWP model information in an objective analysis, which could then be used in the snowfall accumulation algorithm. However, it has long been known that snowfall is significantly underestimated by measurement instruments (Sevruk 1982; Sumner 1988). Indeed, in the context of this study, experiments with objective analyses of precipitation gave snow depths with biases larger in absolute value than those obtained with precipitation forecasts. Although methods of correcting the measurements have appeared in the literature, they vary with gauge type and, at present, no information on gauge type is included in the real-time reports used here.
The snow depth analysis can be simply converted to snow cover by choosing an appropriate threshold. Using a 1-cm threshold, the snow depth analysis was verified against SSM/I-derived snow cover, and the results appear in section 4b. The interpretation of the results was complicated by the apparent errors in the SSM/I snow. Nevertheless, some conclusions may be drawn. In the Northern Hemisphere, the analysis is able to detect all but 6% of the snow cover detected by SSM/I. In western North America, where snow cover is frequently confined to mountain ranges or plateaus, the analysis qualitatively agrees with SSM/I snow. However, over the mountain ranges of Europe and western Asia, the analysis frequently failed to detect snow cover. This error appears to be due to temperatures that are too warm—the result of an inappropriate specification of soil temperature that affects primarily high elevations in the NWP model. Analysis false alarm ratios were generally high, but these could usually be traced to problems with SSM/I snow over densely forested areas, near coasts, where the snow cover is thin (i.e., less than 5 cm), and where the snow is melting. It should be pointed out, however, that the SSM/I algorithm is generally good in mountainous areas where snow cover is typically undetected by the synoptic station network. It would likely benefit the analysis to incorporate SSM/I snow in these areas.
The snow depth analysis was also verified against the NOAA visible product for a 31-month period (section 4c). The results are more difficult to interpret due to the asynoptic nature of the visible product and the large difference in spatial resolution between the visible product in digital form and the analysis. Nevertheless, this part of the study served to identify a problem with excessive snow over the Tibetan Plateau. The problem is due to a combination of excessive model precipitation during the winter months, below-freezing temperatures, and a lack of real-time in situ measurements of snow depth from the region. With the Tibetan Plateau excluded, an overall correlation of 0.86 was obtained between analysis and observation, with the correlation for individual weeks occasionally exceeding 0.91. It may be possible to obtain better agreement if a higher-resolution visible product were used, such as the 24-km resolution snow cover product currently being tested by NOAA/National Environmental Satellite, Data and Information Service.
Apart from its use in NWP forecasts, the snow depth analysis scheme may be useful in some other fields of endeavor. Since snow cover is an important climate variable, this product may be useful to identify and monitor climate change. Developers of snow water equivalent or snow depth algorithms for passive microwave satellite data may be able to use the snow depth and density information to calibrate their algorithms. In the context of NWP, the snow depth and density may be useful to verify winter precipitation forecasts. Undoubtedly, the utility of the product could be greatly enhanced if precipitation were more accurately known. Future work will be directed at improving the quality of precipitation information, either by improving the forecasts or by developing an accurate analysis of precipitation, and at introducing data from passive microwave sensors.
I am grateful to R. D. Brown for many stimulating discussions on snowpack modeling and for suggestions on an earlier draft of this paper. Special thanks to H. L. Mitchell for his encouragement and for his careful review of the manuscript. I thank the three anonymous reviewers for helpful information and suggestions, R. Ferraro for his help validating the SSM/I algorithm, and N. Wagneur for his help with the SSM/I data.
Corresponding author address: Bruce Brasnett, Canadian Meteorological Centre, 2121 Trans-Canada Highway, Third Floor, Dorval, PQ H9P 1J3, Canada.