Introduction

In attempting to use satellite multi-spectral remote-sensing data to map snow and estimate snow characteristics over remote and inaccessible areas, we. are faced with several problems: (1) We must distinguish snow from other surface cover and from clouds; (2) We must compensate for the effects of the atmosphere and rugged terrain; (3) Our space-borne radiometers typically measure reflectance in a few wavelength bands, but for climate modeling we are interested in snow’s reflectance throughout the solar spectrum. Thus, our objective is to map snow and classify it into albedo categories from an existing satellite, the Landsat Thematic Mapper (TM), We also wish to minimize the use of ancillary data. Although digital elevation data may be used to simulate typical conditions for a satellite image, precise registration of an elevation data set with satellite data is often impossible.

The term “mapping” means distinguishing snow from other surfaces or clouds. “Classification” means approximation of snow grain-size and contaminant amount, such that snow spectral albedo can be calculated with a radiative-transfer model.

Spectral Signature of Snow in Thematic Mapper Bands

Table I, which uses data from Reference Markham and BarkerMarkham and Barker (1986), specifies the wavelength bands and saturation radiances for the Landsat-4 and Landsat-5 Thematic Mappers, For the solar part of the electromagnetic spectrum, the values of the exo-atmospheric solar irradiance at the mean Earth-Sun distance integrated over the wavelength bands are also given. The solar-irradiance data used are from Reference Neckel and LabsNeckel and Labs (1984) and Reference IqbalIqbal (1983). In the last column of the table, the sensor saturation radiance is expressed as a percentage of the exo-atmospheric solar irradiance. If the product of the planetary reflectance and the cosine of the solar zenith angle exceeds this value, the sensor will saturate in this band.

The reflectance of snow can be modeled as a multiple- scattering radiative transfer problem, where the scattering properties of the grains are mimicked by some sort of “equivalent sphere” and near-field effects are assumed unimportant. Reference WarrenWarren (1982) has discussed these issues in his review of the optical properties of snow. In the visible wavelengths (TM bands 1 and 2 especially), ice is highly transparent and snow reflectance is insensitive to grain-size but sensitive to modest amounts of absorbing impurities. In the near infra-red (TM4), ice is moderately absorptive and snow reflectance is insensitive to absorbing impurities but sensitive to grain-size. In TM5, ice is more strongly absorptive than in the shorter wavelengths and is more absorptive than water. Snow reflectance in this band is sensitive to grain-size only for very small radii; for most snow, reflectance will be near zero. Both ice and water clouds, however, will be appreciably brighter than snow, allowing for snow/cloud discrimination. Table II shows reflectance for pure snow in the TM bands, and Table III shows reflectances for ice and water clouds for comparison. These were computed by a delta-Eddington approximation (Reference Nussenzveig and WiscombeWiscombe and Warren 1980), which is one of the class of two-stream equations discussed in a later section. It is particularly useful for forward-scattering media, such as ice grains. The refractive-index data are Reference Hale and QuerryHale and Querry’s (1973) for water and Reference WarrenWarren’s (1984) for ice, and the Miescattering approximation is from Reference Nussenzveig and WiscombeNussenzveig and Wiscombe (1980). When absorbing impurities, dust or soot, are present in the snow, reflectances are reduced most in TM band I, moderately in TM2, and less so in TM3 (Reference Warren and WiscombeWarren and Wiscombe 1980, Reference Grenfell, Perovich and OgrenGrenfell and others 1981). The effect of moderate concentration of impurities in TM bands 4, 5, and 7 is negligible. Unfortunately, when impurities are present, the reflectance in the visible bands is also sensitive to grain-size. We did not calculate Mie-theory adjustments for contamination to the pure snow reflectances, because of the uncertainties about how to treat impurities that are embedded in the ice grains (Reference BohrenBohren 1986). Instead, we note that the measurements of Reference Grenfell, Perovich and OgrenGrenfell and others (1981) show that albedo in the visible-wavelength TM bands is reduced by moderate contamination about 0.05 in TM1, 0.03 in TM2, and 0.02 in TM3.

The snow-albedo model calculates reflectances only for monochromatic wavelengths, and thus cannot be directly used to obtain band-integrated reflectances for the TM, because both the refractive index of ice and the spectral distribution of the incoming irradiance vary over the bands. For the present investigation, however, we assume we can effectively mimic band-averaged values by monochromatic optical properties.

The parameters that depend directly on wavelength and grain-size are the snow’s single-scattering albedo and asymmetry parameter g _s; for finite-depth snow, the extinction efficiency Q _ext is also needed. These dependencies can be approximated with empirical equations as functions of the grain radius r.

Table IV shows the values of the coefficients for the seven reflective TM bands. Any coefficients not significant at the 10⁻⁶ confidence level are set to zero.

Use of Digital Elevation Data in Radiation Calculations

Most radiation calculations over terrain are made with the aid of digital elevation grids, whereby elevation data are represented by a matrix. In the USA, these are available as “Digital Elevation Models” from the US Geological Survey (Reference Elassal and CarusoElassal and Caruso 1983). The 1:250 000 scale 1° × 2° quadrangles for the entire USA are available at 63.5m grid resolution (0.01 in at map scale), and the 1:24 000 scale 7.5 min quadrangles are available at 30 m resolution for parts of the country.

From the DEMs, the slope angle s and exposure azimuth E from south can be calculated. We consider a right-handed coordinate system with x increasing toward south and y increasing toward east. The partial derivatives of elevation in the x and y directions are computed from finite differences, and the equations for S and E are given below. Note that the signs of the numerator and denominator allow E to be uniquely determined over [—π,π].

The most important variable controlling the incident radiation on a slope in mountainous terrain is the local solar illumination angle which is determined from the slope orientation, the solar zenith angle on a horizontal surface θ_s, and the Sun’s azimuth Φ₀.

Reference Dozier, Bruno and DowneyDozier and others (1981) described a fast method to determine which points in a DEM are hidden from the Sun by a local horizon for a given solar azimuth. For each point, we calculate the horizon angle H(_Φ) in any desired direction Φ. The horizon can result either from “self-shadowing” or from adjacent ridges. For an unobstructed horizontal surface H(Φ)= π/2.

For diffuse irradiance, only a part of the overlying hemisphere is visible. The “sky-view factor” V _d is the ratio of the diffuse-sky irradiance at a point to that on an unobstructed horizontal surface, i.e. 0 < V _d ≤ 1. We define an anisotropy factor η _d such that η(θ,Φ)L ^↓ = L(θ,Φ), where L ^↓ is the average down-welling radiance on a horizontal surface. Therefore η _d is normalized such that its hemispheric integral projected on to a horizontal surface is π, so for isotropic diffuse irradiance η _d = 1. V _d on a slope S with exposure E is given by

If diffuse irradiance from the sky is isotropic, the inner integral above can be evaluated analytically. We therefore often use the approximation

By a similar “terrain view factor”, it is also possible to, account for reflected radiation from the surrounding terrain, but the formulation is more complicated because the isotropic assumption is almost always invalid (Reference ArnfieldArnfield 1982). Anisotropy results from differing illumination on the surrounding terrain and from geometric effects between the point and the surrounding terrain, even if the surfaces are Lambertian

η _t, accounts for the anistropy of the reflected or emitted radiation. The limits of integration for the inner integral are from the horizon downward to where a ray is parallel to the slope:

In the up-slope direction, cos(Φ – E) is negative, so Ψ_(Φ) < π/2, In the down-slope direction, cos(Φ – E) is positive, so ψ _(Φ) > π/2. Across the slope, Ψ_(Φ) = π/2. Rigorous calculation of V _t is difficult because it is necessary to consider every terrain facet that is visible from a point in order to calculate η _t. No-one has yet done it. We therefore note that V _d for an infinitely long slope is (1 + cosS)/2 and use the approximation

These terrain specifications allow us to specify boundary conditions for an atmospheric radiation calculation over mountainous terrain. However, any mapping or classification algorithm that requires that satellite data be precisely registered to accurate digital elevation data is severely constrained and probably doomed to failure. Digital elevation data in mountainous areas are often of poor quality, with considerable noise from the digitization process, and the differencing operations needed to calculate slope and exposure amplify this noise (See Figs 1 and 2). Thus DEMs can be used to simulate the effects of the combination of atmosphere, terrain, and surface reflectance, and they can be used for radiation models, where the results are integrated over a drainage basin. Moreover, registration between a satellite image and a DEM can be achieved closely enough so that the elevation of a pixel in the image can be known, but an algorithm that needs the slope and azimuth of a given pixel in order to interpret its multi-spectral radiometric signal imposes an impossible requirement, given the poor quality of available DEM data.

Fig. 1. Shaded relief image of the Mount Tom area in the southern Sierra Nevada made from a DEM. The striping in the shaded relief image is from noise created when the DEM is made. Compare this image with Fig.2.

Fig. 2. Thematic Mapper image registered to the DEM in Fig.1. Blue = TM2, green = TM3, and red = TM4. The topographic details are much sharper than in the image made from the DEM in Fig.1. The dimensions of this scene are 15 km × 15 km.

Radiance Above the Atmosphere

The satellite measures upwelling radiance L ^↑ above the Earth’s atmosphere, i.e. at optical depth τ = 0. Interpretation of this signal depends on the interaction between the surface reflectance properties, the terrain, and the atmosphere. Because of the difficulties in dealing with the mountainous terrain, it makes little sense to use a sophisticated atmospheric model. Instead, we use a simplified “two-stream” model (Meador and Weaver I980), which is defined by the following pair of differential equations for a homogeneous plane-parallel atmospheric layer:

where πS ₀ is the exo-atmospheric solar irradiance incident at angle θ₀; μ₀ = cosθ₀; and _a is the atmosphere’s single-scattering albedo. Since the atmospheric layer is homogeneous, the γ values are by definition independent of optical depth. The γ values depend on the single-scattering albedo _a, the scattering-asymmetry parameter g _a, and μ₀. How they are calculated is determined by the particular approximation to the scattering-phase function and the radiation-intensity distribution. Reference Meador and WeaverMeador and Weaver (1980) gave expressions for seven different two-stream approximations.

For a single-layer atmosphere, the solution to Equation (9) depends on the boundary conditions. The usual top- boundary condition is that there is no diffuse irradiance at the top, i.e. L ^↓ ₍₀₎ = 0. At the bottom (optical depth τ₀), however, the terrain effects must be included in the boundary condition, and no-one has yet found a rigorous solution to this problem. Here we resort to an ad hoc method. We assume that the atmospheric transmission and scattering properties are the same as over a flat surface with the same mean regional albedo p _a, so the lower-boundary condition for Equations (9a) and (9b) is

The solutions for the upward radiance at the top L ^↑ ₍₀₎ and the downward radiance at the bottom L ^↓ _(τ0) are:

These equations define reflectance and transmittance for the layer. The other symbols are:

Because of scattering, the relationship between L ^↑ ₍₀₎ and L _(τ0) depends on ρ_g. However, for a thin atmosphere, one can assume it is linear, and we find the relationship by evaluating Equations (11a) and (11b) for ρ_g = 0 and ρ_g = 1.

Now, to find L ^↑ _(τ0) for a surface in rugged terrain, we assume that the same transmission relation holds. The upwelling radiance for a surface with direct albedo ρ_s(μ_s) at illumination angle θ_s, diffuse albedo ρ_d, and albedo of the “surrounding” terrain ρ_g, is the sum of several terms: reflected direct irradiance, reflected diffuse irradiance, and irradiance reflected from adjacent slopes toward the point. The sum of these must be multiplied by cosS to find the mean upwelling radiance projected on a horizontal plane. Therefore, using the simplification in Equation (8) and letting μ_s = cosΦ_s, we find the necessary value to insert L ^↑(τ₀) in Equation (12) above

L ^↑(τ₀) given by Equation (11b). For simplification, we assume that ρ _g = p _s(µ0). This is reasonable for thin atmospheres, because direct irradiance is the largest component of the illumination.

Analysis of Thematic Mapper Signals from Snow

The atmosphere/terrain radiation model described in the previous sections is combined with calculations of the spectral reflectance of snow to simulate radiance at the top of the atmosphere for a range of snow grain-sizes, contamination amounts, and terrain conditions. The atmosphere used in the simulations is the US Standard with a rural background aerosol, 23 km surface visibility, and 50% relative humidity. The optical depths are adjusted to a surface pressure of 650–700 mbar. Table V lists the atmospheric optical properties in the reflective TM bands.

Figs 3 and 4 show a TM image of the southern Sierra Nevada on 10 December 1982. The solar-zenith angle θ₀ is 64.6° and the azimuth Φ₀ is 31.9°. For these values we calculated top-of-atmosphere radiance for 500 randomly chosen locations in the digital elevation grid, covering the range of terrain characteristics listed in Table VI, each for a random grain radius between 50 and 1500 µm, For bands TM1, TM2, and TM3, the simulations included clean, moderately contaminated, and dirty snow.

Fig. 3. Thematic Mapper image of the southern Sierra Nevada. Blue = TM2, green = TM4, and red = TM3. The dimensions of the scene are 85 km × 80 km. It is sub-sampled such that only every sixth pixel is printed, because the size of our image display screen is 512 × 512.

Fig. 4. Same scene as in Fig.3, but blue = TM5, green = TM2, and red = TM4.

The top row in Fig.5 shows theoretical cluster plots of the radiance values from this simulation, paired by some combinations of TM bands. The axes are scaled to correspond to the band-saturation values listed in Table I, From these graphs, it is apparent that the top-of-atmosphere radiance over a snow layer falls into some convenient envelopes for a range of grain-sizes and contamination amounts. The most useful clusters are TM1 and TM4 or TM5 for the shadowed areas, where a threshold brightness distinguishes snow in shadow from other surfaces, and TM2 and TM5, where a threshold brightness distinguishes clouds, vegetation, soil, or rocks from snow in the sunlit areas. Note that most of the values in TM1 saturate, so only a small part of the cluster is on the graph. TM2 and TM3 are redundant. The bottom row in Fig.2 shows actual cluster plots from the TM scene in Figs 4 and 5. These can be used to map snow automatically from the corresponding TM image, without registering the image to a digital elevation grid, A threshold value is chosen from the image for the lowest TM1 radiance for snow in the shadows (∼70 W m⁻² μm⁻¹ sr⁻¹ for this midwinter scene). Then TM2 and TM5 are used to distinguish snow from clouds or other surfaces in the sunlight. In this midwinter scene snow is characterized by:

Fig. 5. (Top) Theoretical cluster plots of pair-wise combinations of Thematic Mapper top-of-atmosphere radiances for snow surface corresponding to the scene and the topography in Figs 3 and 4. The axes limits are the saturation radiances for the sensor given in Table I. (Bottom) Actual cluster plots of pair-wise combinations of Thematic Mapper top-of-atmosphere radiances for the southern Sierra Nevada, corresponding to the scene in Figs 3 and 4.

Fig.6 shows a map of the snow-covered area made by these criteria. The colours in the snow-covered area are stretched, so that the brightest whites represent new, finegrained snow, the brownish tints represent older metamorphosed snow, and the purple areas represent vegetation growing above the snow cover.

Fig. 6. Automatic mapping of snow-covered area in the southern Sierra Nevada scene using the clusters in Fig.5.

Because of the low saturation value in band TM1, it is not feasible to determine snow-contamination amount. The differences between dirty snow and clean snow are not great enough to overcome the other sources of signal variation in mountainous terrain.

Conclusions

The multi-spectral signature from the Landsat Thematic Mapper can be used effectively to map snow and classify it into albedo categories without requiring that the satellite image be registered to digital elevation data. Thus it is possible to use such data for snow-surface energy-balance models in alpine areas.