Introduction
Single-look microwave satellite observations of snow-covered watersheds are being used to produce estimates of snowpack water equivalent (e.g. Reference Foster, Rango, Hall, Chang, Allison and DiesenFoster and others, 1980; Reference Kunzi, Patil and RottKunzi and others, 1982; Reference GoodisonGoodison, 1990; Reference Hallikainen and JolmaHallikainen and Jolma, 1992). However, grain-size, snowmelt and atmospheric effects can shift the geophysical algorithm and scatter the estimates. Multi-temporal looks over the winter season were examined to see if snow-resource estimates can be improved using a sequence of satellite observations. The observations were interpreted using a multiple-scattering model, to assess the model's utility in helping to correct water-equivalent algorithms for the effects of physical factors.
Fig. 1. Target areas in the Canadian plains.
Satellite Observations
Nimbus 7 satellite observations were analyzed for eight winters, along with ground-based hydrometeorological observations. The sample was split, with the first four winters used for model calibration and the last four winters reserved for model testing. To date only the winter 1984—85 has been tested. The apparent brightness temperature was computed for each band, each brightness an average of the pixel values weighted by included target area. The target areas are located in the Canadian plains (Fig.1).
The development of snow water equivalent over time was simulated from surface observations of daily snowfall and maximum air temperature (Reference WankiewiczWankiewicz, 1992). Daily observations were averaged from three to four climatological stations near each watershed. Hourly air-temperature observations were obtainedfrom the closest synoptic station to each basin. Daily soil temperature (2 cm depth) was obtained from the closest agricultural field station. The hourly data were used to categorize the snow state at the time of each satellite orbit. As an example, microwave spectral differences are categorized in Figure 2. In the diagram a transition state between warm and cold allows for uncertainty in measurement,spatial variability or residual wetness; indicated either by a current air temperature Ta of between -1°C and + 1°C; or by positive Ta on the previous day. Although cold snow is indicated by current Ta —1°C or previous day Ta < 0°C, it could include instances of refrozen snow from earlier melt events.
Methods for Satellite Snowpack Monitoring
The conventional method for satellite snowpack monitoring is to develop microwave brightness algorithms, for
Fig. 2. Example of snow categorization of microwave observations.
Fig. 3. Linear algorithm for snow water equivalent from Nimbus 7 SMMR observations.
snow water-equivalent. As example of this approach, a linear regression for Antler River target area is shown in Figure 3 for the calibration period 1978-79 to 1981-82.
A more physically based method is to use a realistic model of microwave-snow interaction. Analogous to modelling watershed run-off, modelling brightness involves setting parameters to either physical inputs, if known; or making simulations match calibration period observations. The multiple-scattering model RADFIT was coded to output time series of brightness temperatures from input time series of hydrometeorological variables. The Appendix describes RADFIT, designed to run on a PC486 microcomputer. It utilizes matrix doubling to compute volume scattering and transmission and simple rough boundary routines to compute surface scattering. Currently, the model is coded for sparse-media Rayleigh scattering. However, the radiative-transfer formulation used in the model allows straightforward substitution of the dense-media low-frequency scattering expressions given in Reference Wen, Tsang, Winnebrenner and IshimaruWen and others (1990).
Model Calibration Using Binned Values of Snowpack Mass
RADFIT is a forward algorith, but can be iteratively used for the reverse problem of deriving snowpack parameters from microwave observations. This was done by first binning the data into five water-equivalent cases, each case representing over 30 observations. Also, daily hydrometeorological variableswere assembled, period averages of which are shown for each basin in the upper section of Table 1 (the middle section shows the corresponding period-average brightnesses). Then, for each model iteration, one or more snowpack or soil parameter was changed, to produce a closer match between computed and observed brightnesses. In practice, a fast version of the model was initially iterated, with eight-angle discretization over the polar angle and a uniform snow temperature. Final tuning of the iterations was performed with 14 polar angles. Error in using a single layer was assessed by comparing the results with those from a ten snow-layer model.
Final values for fitted snow and soil parameters are shown in the last section of Table 1, for all three basins. The derived grain-size does not change much between target areas, in spite of the great sensitivity of brightness temperature to its magnitude. The derived dirt factor, an effect of snow impurities on microwave absorption, shows a greater variability. As shown in the Appendix, snow transmission is temperature dependent, with a sensitivity assumed to vary with contaminant concentration. Direct measurements of impure ice absorption (e.g. Reference Matzler and WegmüllerMätzler and Wegmüller, 1987) are being reviewed for possible incorporation in the next model version. For soil
Table 1. Parameters for the multiple scattering model (average period values are given for the variables)
parameters, model fitting produced small soil roughness with considerable polarized reflectivity mixing (defined in Appendix), indicating a smooth surface with small polarization. The fitted soil matrixproperties indicate soil saturated with ice/water; with an unfrozen fraction of 10—15% in spite of the sub-zero soil temperatures shown in Table 1.
All four frequency-polarization combinations are solved simultaneously. The closeness of the resultant fit is shown for Antler River in Figure 4, a result typical of all three basins (note that each basin has been fitted with its own parameters listed in Table 1). The figure shows that when a close fitis made with the current version of the model to the rapidly varying 37 GHz channel (primarily by adjusting grain and dirt factors), a flat response results at the lower frequency. The slow increase in brightness at 18 GHz is a result of the derived radiometrically cold soil. The two channels cannot beindependently fitted because scattering and absorption coefficients are linked via the scattering andice absorption-frequency dependence. Nor can exact simultaneous fits be made at the y axis (byadjusting soil parameters) because boundary reflectances are linked via the Fresnel equation and water-relaxation frequency dependence.
Simulation of Brightness Temperature Snowpack Variations
The variation of brightness temperature as a function of snowpack was investigated for the 1984—85 winter season, the first and only year simulated so far. Figure 5 shows the observed changes in Antler River microwave spectral difference as a time series represented as solid circles connected in a sequence moving (generally) from left to right. The observations are compared with model results for bothlong- and short-term brightness changes.
For interannual comparison, the linear algorithm from Figure 3 (1978-82) is duplicated in Figure 5 asthe dashed line. It lies below the 1984-85 observation series. Testing of additional years may revealwhether this shift is a temporal fluctuation or a long-term trend. At any rate, accuracy for late season forecasting in this particular year would improve if the algorithm was translated upwards to match early-season multi-temporal observations.
The complete set of input variables was used in the multiple-scattering model for winter 1984—85, brightness temperatures being computed for each of the 33 available satellite orbits. An interannual shift is also apparent in the simulated time series as well (open diamonds in Figure 5). In addition, its seasonal trend differs in slope from the other two time trends, because of the effect on brightnessdifference of an optimizing bias in favor of 37 GHz, described above. The use of more sophisticated routines is currently under investigation (e.g. dense-media theory, impure snow-ice absorption observations, multi-layer atmosphere) to see if a better fit can be produced for the two-frequency set.
Both the observations and simulation time series show short-term fluctuations unrelated to snowpack mass. Overall, the fluctuations in the modelled series are also seen in the observation series, although not always of the same magnitude, indicating that the model subroutines
Fig. 4. Fit of multiple-scattering model results to observations.
do represent real physical processes. A detailed examination of both series, however, reveals three model points whose behaviour is not reflected by that of the corresponding observation points. In these cases it was found that the snowpack had been previously exposed to either a melt cycle or a rainfall event, and that the observed brightness was lower than expected. Grain and ice-layer development is not represented in the model.
Fig. 5. Simulation of seasonal brightness temperature changes.
Conclusions
The following conclusions on satellite procedures for snowpack monitoring are suggested by inspectionof Figure 5.
-
Both linear algorithm and physical snowpack models will have improved accuracy if multi-temporal observations of microwave brightness are utilized. By combining these observations with recent snowpack surface measurements, compensation can be made for errors from inter-year changes in previously fitted parameters.
-
The multiple-scattering model appears to reproduce short-term fluctuations in microwave satellite brightness temperature, indicating that microwave interactions are properly modeled at least on this time-scale. Further testing is required to see if the corrections for these effects can be inverted, to improve snowpack measurements made in the presence of these fluctuations.
Appendix
Multiple-Scattering Model Radfit
The scattering model comprised three modules: parameter input, computation and output display. The input module provides a selection of interactive and file formats for iterative testing of one or more input variables. The computation module provides multiple loops through the model kernel, computing brightness for a set of observation cases and frequencies. The output module allows visual analysis ofdifferent microwave components, or quick comparison between observation and simulation.
The model kernel is computed in the sequence: electromagnetic properties, multiple volume scattering,rough boundary scattering, layer-scattering operators, emission sources, and terrain brightness from space. The electromagnetic properties are computed from physical models of the media. Dielectric constants are derived for soil water from the Debye equation, soil layer from the mixture of Reference Dobson, Ulaby, Hallikainen and El–RayesDobson and others (1985; extended here to include an ice phase) and dry-snow real part from observation (Reference Ulaby, Moore and FungUlaby and others, 1986, fig. E26). Snow absorption and scattering was calculated from sparse-media Rayleigh scattering. The dielectric constant of the ice component is esi = 3.15 -j × esi″, where esi″ is a function of both frequency f and snow temperature Ts. For given Ts, esi″ = dirt × (C/f + D × f^e); where dirt is an ice-impurity factor and C, D and e are obtained from Reference Ulaby, Moore and FungUlaby and others (1986, fig. E.3 for Ts = -1°C and -20°C and for dirt = 1). This form of f dependence is identical to that of Reference Matzler and WegmüllerMätzler and Wegmüller (1987). esi″ for other Ts are obtained using the interpolation ((A + B) × exp(-kTs) + B), by fitting A and Β to known values and k to Ts functions in Reference Ulaby, Moore and FungUlaby andothers (1986, fig. E4).
Snow-layer temperatures are linear interpolations of values at three levels: surface, base and mid-snow. The first two are set equal to air and soil temperature, Ta and Tg; the latter is derived from a finite-difference conduction routine from the Ta and Tg time series. Atmospheric optical thickness τ is computed for water vapor using a single equivalent layer (Gloersen and Barath, 1984); and for oxygen using observations in Reference CraneCrane (1981; shown in Reference Ulaby, Moore and FungUlaby and others, 1981, fig. 5.7b). Precipitable water was derived from the regression with dew point given in Reference ReitanReitan (1963).
Multiple volume scattering is computed by integration over m refracted angles in the snow, i.e. m - 1 equally spaced cosines and one cosine from the view angle θi. In each of k equal water-equivalent layers, the transmission and reflection (TR) operators are derived from the matrix-doubling operator (Reference Grant and HuntGrant and Hunt, 1969). Multi-layer TR operators are separately derived for energy flow in both the upward and downward directions, by reversing matrix multiplication order. Layer effects on volume scattering can be computed from profiles of absorption, scattering albedo and temperature. For simplicity, inter-layer boundary scattering is ignored, except at the air/snow and snow/ground interfaces.
Snow-surface scattering is computed by a five-fold Gaussian-Legendre quadrature at each of the m angles. This way, the correct angular weighting is given to each reflectivity, especially near thecritical angle. Flat-boundary reflectivities gv and gh (vertical and horizontal polarization) and the refraction angles, come from Fesnel's and Snell's equations. Rough boundary reflectivities Gv and Gh come from Reference Wang and ChoudhuryWang and Choudhury (1981):
Q is a polarized reflectivity mixing parameter, h is a boundary roughness and θ is the incident angle at the upper boundary.
Surface and volume scattering are derived from ray tracing (Reference LeaderLeader, 1975). The phase matrices for surface volume scattering are computed for the two surfaces, in both upwards and downwards directions. Multiple boundary reflections are then included in scattering operators for each of four sources, in both directions. The emission sources are: upwards ground emission, from Ts observedat 2 cm depth; downwards emission from the sky “layer” computed by
from the effective atmospheric temperature Te and cosmic source Textra, and both upwards and downwards snow emission. Snow emission is computed either from a memory-extensive layered-snow routine (Reference TwomeyTwomey, 1979; requiring numerous layers for accurate results) or from Kirchhoffs law for a single-snow layer.
For the model runs described in this paper, a relatively crude single-layer routine was used to compute sky radiation. The effective temperature Te was arbitrarily set to (Ta — T0), with T0 = 2 K. Since then, a more accurate 64-layer routine has been programmed; consequently the error in the former approach can be assessed. The multi-layered atmospheric routine is valid to 300 GHz, and follows Reference RosenkranzRosenkranz (1975; in Reference Ulaby, Moore and FungUlaby and others, 1981, p. 276–78) and Reference WatersWaters (1976; in Reference Ulaby, Moore and FungUlaby and others, 1981, p. 271–74). Both routines were compared for January pressure, temperature and moisture profiles given in Reference SissenwineSissenwine (1969) at both 45° and 60°N; over a bare ground surface of temperature Ta and emissivity 0.77. It was found that, compared to the multi-layer model, the single-layer model overestimated brightness temperatures from space: by 0.9 Κ and by 0.5 K, at 45° and 60°N, respectively, for 18 GHz; and by 3.1 Κ and 1.1 K, respectively, for 37 GHz.
The resulting terrain brightness at satellite altitudes follows after reduction by the loss factor = exp(τ × secθi), and augmentation by upwards atmospheric emission Te × (1 — 1/L). Finally, the brightness TAP at each polarization is computed by averaging snow terrain and bare-ground values (TAPSNOW and TAPGND):
where cover is the fractional snow cover for the pixel. The kernel is re-computed for each new frequency and case combination.
