Introduction

Satellite passive-microwave observations of snowpacks have the capability of estimating parameter characteristics of the snowpack. Reference EnglandEngland (1974, Reference England1975), Reference Chang., Gloersen, Schmugge, Wilheit and Zwally.Chang and others (1976) and Reference Tsang. and Kong.Tsang and Kong (1977) used the radiative-transfer theory (Reference ChandrasekharChandrasekhar, 1960) to model the electromagnetic response of snow. Reference Chang., Gloersen, Schmugge, Wilheit and Zwally.Chang and others (1976) developed a scattering model based on Mie scattering theory that predicts the brightness-temperature response to variations in snowpack thickness and grain-size. The intensity of the microwave radiation emitted from a snowpack depends on its physical temperature, grain-size, density and the underlying soil conditions. Snow grains scatter the upwelling emission and thus the amount of scattering will increase with snow depth, resulting in a decrease of brightness temperature measured by the microwave radiometer. Extensive experimental work to measure the interaction of microwave radiation with the snowpack has been undertaken with a passive-microwave instrument mounted on a tower 15 m above ground (Reference MätzlerMätzler, 1987). The radiometer measures effective physical temperatures at 4.9, 10.4, 21.35 and 94 GHz, and estimates the corresponding emissivities. While these measurements have clearly shown the dependence of brightness temperature on the snowpack properties, it is the spatial variation of the snowpack within the SSM/I footprint that is not a priori known and that is difficult to estimate. Other studies have been conducted with in-situ snowpack observations to derive algorithms to determine snow water equivalent from satellite passive-microwave observations (Reference Hallikainen and Jolma.Hallikainen and Jolma. 1986; Reference McFarland, Wilke and Harder.McFarland and others, 1987; Reference Hall.Hall and others, 1991; Reference Goodison., Walker., Choudhury, Kerr, Njoku and Pampaloni.Goodison and Walker, 1995; Reference Rott, Nagler., Choudhury, Kerr, Njoku and PampaloniRott and Nagler, 1995; Reference Josberger, Gloersen, Chang and Rango.Josberger and others, 1996). The fraction of the sensor footprint which is covered with forest adds another complication to the development of a suitable algorithm. The vegetation canopy can be represented as a dielectric mixture consisting of dielectric elements (leaves, stalks and branches) embedded in a matrix of air (Reference Ulaby and Jedlicka.Ulaby and Jedlicka, 1984). A simple model to separate the effect of the forest from the effect of snow was proposed for the SMMR (Scanning Multifrequency Microwave Radiometer) 37 GHz channel by Reference Hall, Foster and Chang.Hall and others (1982) and an algorithm to correct for the absorption of the snow signal by the forest cover was subsequently proposed by Reference Chang and Chiu.Chang and Chiu (1991).

In the development of snowpack algorithms, the snowpack thickness or water equivalent within the satellite footprint is generally derived from point measurements. Because of the size of the satellite footprint, which is approximately 30 km for the 37 GHz channel, even a large number of point measurements on the ground may not correctly reflect the snowpack properties across the entire pixel, typically, there are only a few point measurements within the satellite footprint (Reference Künzi, Palil and Rott.Künzi and others, 1982; Reference Josberger, Ling, Campbell, Gloersen, Chang and Rango.Josberger and others, 1989; Reference Hall.Hall and others, 1991; Reference Hallikainen. and Jolma.Hallikainen and Jolma, 1992; Reference Goodison., Walker., Choudhury, Kerr, Njoku and Pampaloni.Goodison and Walker, 1995). While the resulting correlation between the SWE and the satellite observations can be high, the snowpack can be very variable across a 30 km footprint. Algorithms developed from snowpack information averaged over spatial scales comparable to the satellite observations would be more reliable and could validate other algorithms.

The airborne snow water-equivalent measurements carried out by NOHRSC provide an alternative data source that is nearly on the scale of the satellite observations. Throughout the winter, NOHRSC derives the average snowpack water equivalent (SWE) along flight lines that are approximately 10 km long, by measuring the attenuation of naturally emitted terrestrial gamma radiation by the snowpack (Reference Carroll and Voss.Carroll and Voss, 1984; Reference Carroll and Carroll.Carroll and Carroll, 1989). The error associated with these SWE determinations over agricultural vegetation is 8.2 mm of SWE and over forest land use the error is 23.2 mm of SWE. A large number of these flight lines in a region a few hundred kilometers wide, flown in a short period of time, would provide a unique dataset to compare with satellite observations. Such was the case in 1989, when NOHRSC flew 92 flight lines in the Red River Valley region in mid-February and early-April 1989. The region spans the boundary between Minnesota and the Dakotas, has very little relief and the eastern part is forested while the remainder is open. Hence, the airborne measurements combine with the physiographic characteristics of the region to provide a unique opportunity of studying a prairie-type snowpack as categorized by Reference Sturm, Holmgren and Liston.Sturm and others (1995).

The Datasets

The datasets used in this study consist of the NOHRSC flight-line SWE determinations, passive-microwave observations from the SSM/I carried by the F8 satellite of the DMSP and land-use classification as determined from observations by the AVHRR. For the Red River Valley region, NOHRSC flew 92 flight lines in 3 days, 15×17 February 1989. Figure 1 shows the location of the centers of these lines in the study area and the snow water-equivalent map obtained with the gamma-ray instrument. The lower left corner of the area shown is located at 45°N and 99°W. A similar airborne survey was carried out on 1–2 April 1989. However, the April observations were not used in this study, because the surface temperatures as measured by the National Weather Service (NWS) stations in the area, were near 0°C. This resulted in liquid water in the snowpack and little correlation between the satellite observations and the SWE.

Fig. 1. Snow water equivalent map (cm) derived from NOHRSC flight lines that took place 13–17 February 1989 in the Red River basin. Each cross represents the location of the center of each flight line. The contour lines and gray tones show the gridded SWE values derived by kriging the flight line data.

The SSM/I is a seven-channel passive-microwave radiometer with vertically and horizontally polarized channels at 19.35,37 and 85.5 GHz, and a vertically polarized channel at 22.235 GHz. This study used horizontally (H) polarized data, because there is an inverse relationship between snow depth or snow water equivalent and the 37_H GHz brightness temperature as measured by the satellite sensor under drysnow conditions. The 19_H GHz channel can be used to eliminate partly the effects of the snow and ground temperature and the atmospheric quantities such as integrated water vapor and clouds on changes in brightness temperatures. The microwave data were obtained from the National Oceanic and Atmospheric Administration (NOAA) National Environmental Satellite, Data, and Information Service (NESDIS) in latitude by longitude pixels. The swath from 17 February 1989 covered the region of interest; the early morning pass was used to minimize the influence of the presence of liquid water in the snowpack. NWS surface-temperature observations from the region show that the maximum temperatures were approximately –10°C or less and minimum temperatures were approximately –25°c. In order to correct for the effect of the vegetation on the microwave signal, no-snow microwave conditions for the region were obtained from a satellite pass on 14 November 1988. NWS records from the region show no snow on the ground at the time and the snowpack began to accumulate shortly after this date. Hence, these microwave observations reflect the soil conditions immediately before being covered by snow, when maximum and minimum air temperatures were approximately –1°C and –10°C respectively.

The land-use data were obtained from the Conterminous U.S. Land Cover Characteristics Data Set, 1990 Proto-type CD-ROM, as derived from AVHRR observations. The CD is available from the U.S. Geological Survey, EROS Data Center, Sioux Falls, SD. The land-use classification is provided in 157 categories for 1 × 1 km pixels in a Lambert Azimuthal Equal Area (LAZE) projection. Because the forest in the eastern part of the region consists primarily of coniferous trees, which can have a strong effect on the observed brightness temperatures, we simplified the 1 km land use into a binary classification: 0 for no forest and 1 for forest-covered. This CD also gave the elevation for each 1 km pixel.

Dataset Registration

To produce co-registered data at the same spatial scale, we first transformed the satellite microwave data and the flight-line data into the same projection as that used for the AVHRR data, the LAZE projection. Then, a kriging technique was employed to generate 20 × 20 km grids of snow water-equivalent values and corresponding satellite brightness temperatures. Our analysis only used the data within the irregularly shaped area in Figure 1, where most of the flight lines occurred. The same procedure was carried out for the satellite observations obtained on 14 November 1988, the snow-free situation.

For the land-use classification data, which were already in the LAZE projection, we simply averaged all of the 1 km values within each 20 km gridcell. This produced a forest-cover index that varied from 0 to 1, with 0 being completely open and 1 being completely forest-covered. Most of the region is open, consisting of cropland or uncultivated open land, but the eastern part contains significant coniferous forest cover. Similarly, we defined the ruggedness for each gridcell as the standard deviation of the 400 1 km elevation values within each cell. This procedure yielded a total of 164 points that were associated with snow water-equivalent values.

Analyses

We carried out four types of analysis to determine their capability of predicting snow water equivalent from the microwave measurements. These analyses consist of a multiple linear regression, a parameter evaluation for an algorithm containing a vegetation correction, a neural network approach and a general linear model using standard least squares methods. To test the different methods, we included eight variables: the 19_H, 37_H and 85_H GHz channels, the 22_v GHz vertically polarized channel, the difference between the 19_H and 37_H GHz channels, the land-classification values, the ruggedness and the ruggedness squared. Then, we examined closely the relationships between 19_H and 37_H GHz horizontally polarized data and SWE and the forest cover. Figure 2 shows a scatter diagram of the values of 19_H minus37_H GHz channels for both the February snow-covered case and the November snow-free case as a function of the snow water-equivalent measurements derived from the gamma-ray techniques.

Fig. 2. Scatter diagram oj the differences in brightness temperature between the 19_H and 37_H GHz channels as a function of the snow water equivalent from the gamma ray instrument for February 17, 1989, when there was snow and November 14, 1988, when there was no snow.

1. Multiple Linear Regression

A multiple linear regression of SWE on all eight of the variables that combined the snow-free and snow-covered observations, attained an R ² of 0.89. Inspection of the T-ratios for each variable showed that the 19_H and 37_H GHz frequencies were by far the most significant variables in the regression. The remaining variables, including the forest cover, were found not to affect significantly the regression. This curious result will be discussed later. Following Chang and others (1987), who used the brightness-temperature differences between the 19_H and 37_H GHz channels to estimate SWE, a linear regression of SWE on the brightness-temperature differences was performed and yielded the following regression equation

(1)

This linear regression had the same R ² value as the value found in the standard multiple regression using all of the variables and the value of the first coefficient was not statistically different from zero. Equation (1) is very close to the algorithm derived by Reference Chang, Foster and Hall.Chang and others (1987), who found a proportionality constant of 0.48 cm K⁻¹ for snow with a density of 300 kg m⁻³.

2. Parameter Evaluation with a Forest-Cover Correction

Reference Chang and Chiu.Chang and Chiu (1991) proposed an algorithm that incorporated a vegetation correction. The correction is based on the fraction of the pixel that is forest-covered and the 19_H and 37_H GHz brightness-temperature difference as observed during snow-free conditions:

(2)

where TB is the brightness temperature is the brightness temperature for the corresponding snow-free pixel at 19_H or 37_H, and F is the fraction of the pixel cov ered by forest. For each pixel, we computed A using its brightness temperatures in the snow-free and snow-covered ease (TB) and its forest cover fraction (F). To determine the effects of the forest cover, we also computed another set of A values with F set to zero for all cases and then compared the two cases.

Figure 3 shows the results for the two cases; the individual values of A are plotted as a function of the forest fraction for each pixel. The small circles (o) are values with no correction and the pluses (+) are the values with the forest correction. Because of the large number of pixels with no forest, there is little difference between the average value for each case. However, segregating the A values with respect to forest fraction shows the non-negligible effect of the forest cover (plotted along the abscissa in Figure 3). The bold diamond symbol and the bold plus symbol show the average value of A, respectively for the cases without and with forest correction. The corresponding bars give the standard deviation. For F = 0, all values were used in the average for F = 0.1 only those values from pixels with a forest cover of 0.1 or greater were used, etc. With no forest correction, the average value of A increases as the forest fraction increases. However, with the forest correction, the average value of A remains essentially constant at 0.512 cm K for all forest covered conditions. This value is not statistically different from the value of 0.514 cm K found above, or from the theoretical results of Reference Chang., Gloersen, Schmugge, Wilheit and Zwally.Chang and others (1976).

Fig. 3. Variation of the A coefficient defined in Equation (2) as a Junction of the fraction of forest cover. The values computed with forest correction are represented with circles (o), and the values computed with no correction correction are represented with pluses (+). The bold symbols give the average values for each case and the bars give the standard deviation of the averages. See the text for further explanation.

3. Neural Network

We also carried out a neural-network analysis as described by Reference WassermanWasserman (1989) to determine if an artificial intelligence approach would provide a more accurate algorithm. The initial input data consisted of the eight variables: the brightness temperatures from the four passive-microwave channels (19_H, 22_v, 37_H and 85_H GHz), the 19_H – 37_H GHz difference in brightness temperature, the forest cover, the ruggedness and the ruggedness squared. The architecture of the network can be described as 3–1–1: all eight inputs were read into the input layer; the nodes in the input layer were then fed into a layer of three “hidden nodes” which were in turn fed into a one-node output layer. The activation functions between input and hidden layer, and between the hidden and output layer, were similar to the usual “squashed sigmoid” and are therefore constrained to values between –1 and 1. The networks were implemented in MATLAB (Mathworks, 1994), using the neural-network toolbox, and the MATLAB “tansig” function was used for the sigmoid-activation function.

To train the network, the dataset was randomly partitioned into two parts, two-thirds for training and one-third for testing. The technique described by Reference LeeLee (1995) was used to normalize snow water-equivalent values and each of the eight inputs. A similar process was used for the February dataset containing only snow readings and also for trials limiting the inputs for comparison. The networks were trained for 10 000 epochs on the training set with a learning rate of 0.01. Greater values for the training rate resulted in exponential growth of the error function. In training, the error function dropped rapidly until about epoch 1000, and then dropped slowly from then on. Each training/testing cycle was performed ten times so that the variation caused by different initial weights for back propagation could be examined.

Several tests were performed using two different data-sets: a dataset composed of both the data acquired during the February snow period and the November no-snow period, and a second dataset consisting only of the snow-covered case for February. With these two datasets, we performed four different analyses that are presented in Table 1, where we compare the results obtained with the neural network to the R ² values obtained using a standard multiple linear regression (MLR). For the analysis, we started by using the eight initial inputs, then we eliminated two inputs at a time, starting with the inputs that were found to be less significant in the standard multiple-regression analysis. The first two inputs we eliminated were the 22_v and 85_H GHz channels, which did not significantly decrease the R ² values in the training or the test sets when compared to the values using the eight inputs. Then, the two ruggedness parameters were eliminated arid finally the 19_H and 37_H GHz channels were taken out. When the ruggedness parameters were eliminated, the R ² value for the training set decreased by 0.005 but the R ² values for the test set increased by approximately 0.02 when compared to the preceding test sets. Finally, when the only two input parameters that remained were the brightness-temperature difference between 19_H and 37_H GHz, and the fraction of vegetation cover, both R ² values for training set and the test set decreased significantly when compared to the other values. Note that for the third case with 19_H, 37_H and 19_H–37_H GHz, and the forest fraction, MLR will automatically eliminate the difference input, because it is a linear combination of the previous two inputs.

Table 1. R² values from the neural-network analysis and from the multiple linear-regression analysis (MLR) for various data inputs. The upper part is for the case with both snow and no-snow data, while the lower part is for the case with snow only

For the cases using both the snow-free and snow-covered data (left side of the table), the regression coefficient was quite high for both the neural network and the MLR analysis. For the neural-network analysis, the R ² values for the training data decreased slightly as the number of inputs decreased, and the results were slightly lower for the test data. For the MLR study, the R ² values remained nearly constant as inputs were eliminated, showing that the 19_H and 37_H GHz data were dominant in the regression.

When the dataset with only snow was considered (right side of Table 1), the R ² values decreased significantly for both the neural network and the MLR analyses, from approximately 0.97 to 0.85 or less for the neural network results and from 0.89 to 0.60 or less for the MLR results. This behavior is a result of the limited range of SWE in this study, from 6 and 16 cm, and the non-uniform distribution of data over this range. For this ease, the neural network approach greatly outperformed the MLR as evidenced by R ² values that were 0.2 greater. For both analyses, the correlation decreased as the number of inputs decreased.

4. General Linear Models

The combined snow-covered and snow-free data were also analyzed by traditional least-squares regression models (using second order terms and interactions) to give a standard of comparison for the neural network methods. In the cases where all data were used (snow and no-snow), the results are essentially the same as those obtained by the neural network methods (R ² in the range of 0.98 on the full and training sets, and R ² = 0.95 on the test set). The optimal models make use of the 19_H, 37_H and 85_H GHz inputs and the forest fraction. Analysis of the 19_H and 37_H GHz observations from the snow-covered ease only gave an R ² of 0.52. When the forest cov er was added, the R ² rose to 0.69.

In modeling the data which are restricted to cases with snow on the ground and using the full set of variables, the success of the least-squares models again was comparable to that cited above for the neural network models (with some decline in effectiveness when the 85_H GHz variable was deleted). Figure 4 shows the smoothed variation of snow water equivalent vs the difference in brightness temperature 19_H–37_H GHz for the snow-only data. The smoothed plot is distinctly sigmoidal in shape, which corresponds to usual between layers activation function used in neural network models in this study. To construct models of a more traditional (statistical) type that have increased effectiveness, we would probably have to turn to weighted regression models or projection pursuit regression models.

Fig. 4. Smoothed variation of snow water equivalent as a function of the difference in brightness temperature 19_H and 37_H. This sigmoidal curve was obtained with a general linear model.

Conclusions

Snow water-equivalent estimates acquired in February 1989 from the attenuation by snow of gamma radiation measured from aircraft are compared to brightness temperatures measured with SSM/I during the same period over the Red River basin. The snow water equivalent measured with the aircraft gamma-ray instrument are a unique dataset since the aircraft data are at a spatial scale similar to the SSM/I data. These measurements (shown in Figures 1 and 2) represent a limited range of snow conditions (SWE between 6 and 16 cm) which are a limitation to our study.

For a prairie-type snowpack as found in the Red River-area, conventional algorithms based on the 19_H–37_H GHz difference from both the snow-covered and snow-free cases, can be used to estimate the snowpack water equivalent from satellite passive-microwave observations. An algorithm that includes the forest-cover fraction for each pixel further improves the SWE estimates. A neural-network approach can reproduce the observations to a high degree of accuracy, as reflected by R ¹ values of 0.98. Similar results were found using a general linear model. One has to bear in mind that this experiment, which uses only a limited range of SWE shows great promise for determining accurate SWE estimates from passive-microwave observations. The technique should be further refined by expanding the range of SWE and including grain-size measurements.