General framework

The radiation balance of a snowpack in sloping terrain is given by

where Q* is the net all-wave radiation, Κ ↓ is the total shortwave irradiance, L ↓ is the total longwave irradiance, L ↑ is the snowpack thermal exitance, and α is the albedo of the upper snow surface.

The two-stream radiative transfer model devised by Zdunkowski and others (1982) is used to calculate shortwave and longwave downward fluxes on unobstructed horizontal surfaces. This model handles absorption by H₂O, CO₂ and O₃, aerosol and cloud droplets as well as multiple scattering by air and particles. The O₃ distribution is taken fromReference Craig Craig (1965) for a midlatitude situation, H₂O is obtained from radiosonde ascents, and the CO₂ level is acquired at a nearby meteorological station. Partial multi-layered cloudiness is also treated. The results are then used in equations designed to take into account terrain effects, thus enabling to calculate both shortwave and longwave fluxes on sloping surfaces.

Individual components

In order to calculate incoming fluxes or to determine outgoing fluxes from remotely sensed data in high-relief environments, information characterizing surface topography and atmospheric conditions is a prerequisite. Topographic data can be used effectively either to compute shortwave or longwave irradiance on slopes of varying gradient and aspect, or to correct for atmospheric effects(atmospheric transmission and path radiacne varying exponentially with terrain height (Reference Woodham and LeeWoodham and Lee, 1985)) on remotely sensed imagery when deriving physical quantities such as surface reflectance. Atmospheric conditions, on the other hand, are also important to consider in the derivation of incoming and outgoing fluxes. Meteorological parameters such as air temperature and relative humidity, atmospheric aerosol and liquid water content control to various degrees of importance energy receipts at the surface and the satellite platform. In the proposed model, data input primarily consist of digital terrain data, satellite imagery, and radiosonde data.

Shortwave irradiance

Shortwave irradiance (direct and diffuse) on unobstructed horizontal surfaces is computed by dividing the solar spectrum into four intervals (0.28–1.00 μm, 1.00–1.53 μm, 1.53–2.20 μm and 2.20–6.00 μm) making it possible to handle spectral overlap effects for atmospheric constituents (Reference Brühl and ZdunkowskiBrühl and Zdunkowski, 1983). Direct, diffuse sky and total irradiance values for unobstructed horizontal surfaces are determined for every possible terrain height within a study area through adjustment of exponential functions (Reference Sjoberg and HornSjoberg and Horn, 1983; Wood-ham and Lee, 1985). Reference Duguay and LeDrewDuguay and LeDrew (1992) have shown the appropriateness of exponential functions to estimate these parameters, as well as effects due to path radiance and atmospheric transmission in remotely sensed images of mountainous terrain.

The selection of an algorithm for the computation of shortwave irradiance on surfaces of varying slope and aspect is based on an earlier investigation by Reference IsardIsard (1986) in the east slope of the Colorado Front Range. Following Isard's recommendations, Reference HayHay's (1983) model is used together with results obtained from the two-stream model in order to derive total shortwave irradiance on sloping surfaces. The anisotropy of diffuse sky illumination is considered in this model. The equation can be written as

In the equation, the symbol between square brackets [z] indicates that the value of the variable varies with height. Its value is calculated for all terrain heights within the study area. Briefly, the variable S _h defines the direct irradiance to an unobstructed horizontal surface, computed from the two-stream model. The variable i _h is the angle of incidence to a horizontal surface, and is simply computed as the cosine of the solar zenith angle. Values for variables D _h and K _h are also derived from the two-stream model. The variable k ^′ is the anisotropy index which is utilized to separate diffuse sky irradiance into isotropic and circum-solar components. Values for slope (s), angle of incidence (i) and shadow (φ = 0) are generated using the PRODEM software package (Reference Duguay, Holder, LeDrew, Howarth and DudychaDuguay and others, 1989). Finally, a regional albedo (α) can be specified in the equation for snow-covered surfaces in order to take into account the contribution of terrain-reflected radiation to the total irradiance.

Longwave irradiance

Total longwave irradiance on unobstructed horizontal surfaces is also calculated using the two-stream method of Zdunkowski and others (1982). The calculations are carried out by dividing the spectrum into two broad spectral regions. Longwave fluxes in the first region (3.50–8.75 μm) and 12.25–100.00 μm) are calculated using the emissivity model, while those in the second region (8.75–12.25 μm) are computed with the window model.

In order to determine longwave irradiance (atmospheric emittance) on all surfaces in the study area, a variation of the algorithms presented by Reference Marks and DozierMarks and Dozier (1979) and Olyphant (1986) is utilized. Using results obtained from the two-stream model, the equation can be written as follows:

where L ↓ [z] is the longwave irradiance on sloping terrain at altitude z, L ↓_h is the longwave irradiance on an unobstructed horizontal surface, V _f is the sky view factor, ε_s is the emissivity of the surface, σ is the Stefan–Boltzman constant (5.67 x 10⁻⁸ W m⁻² Κ⁻⁴), and T _s is the temperature of the terrain. Term 1 of the righthand side of the equation represents the portion of the longwave radiation emitted by the atmosphere while the righthand side is the longwave emission from the surrounding terrain. In contrast to models which assume average (regional) surface temperatures and emissivities in the determination of the contribution from the surrounding terrain, the model makes use of the surface brightness temperature measurements obtained from Landsat TM thermal band 6.

For the calculation of longwave irradiance, very accurate determination of surface temperatures and especially emissivities is not always essential since the contribution of the second term in Equation (equ.3) is usually minimal. This is because sky view factors for surfaces are often greater than 0.90, even in high relief environments.

Surface albedo

An estimation of surface albedo from remotely sensed reflectance measurements can be derived when the following three factors are considered (Reference BrestBrest, 1987): (1) the spectral reflectance of the surface of interest; (2) the spectral distribution of the irradiance; and (3) the wavelength region of the discrete spectral bands. In our model, the spectral distribution of solar irradiance is considered to be within the 0.28–6.00 μm range.

The determination of albedo from reflectance measurements requires that the spectral reflectance curve of the material being sensed be divided into segments of uniform reflectance. Each segment of the curve can then be represented by a band measurement located within the range of the segment. In this approach, a weighting scheme is assigned to a particular category based on the spectral reflectance curve characterizing it (see Reference BrestBrest, 1987; Reference Duguay and LeDrewDuguay and LeDrew, 1992, for details). Two weighting schemes are devised for snow-covered surfaces by partitioning a typical reflectance curve into four regions: 0.28–0.725 μm, 0.725–1.00 μm 1.00–1.40 μm and 1.40–6.00 μm. Snow reflectance is very high in the visible portion of the electromagnetic spectrum, but decreases rapidly in the near-infrared and mid-infrared wavelengths. In general, the albedo of snow-covered surfaces can be estimated as

However, when snow saturates in TM bands 1, 2 and 3, TM band 4 is used to estimate reflectance (ρ) for the 0.28–0.725 μm segment. Through the examination of typical snow-reflectance curves (e.g. Reference DozierDozier, 1984;Reference Hall, Chang and Siddalingaiah Hall and others, 1988), the reflectance of snow in TM band 2 was evaluated as 1.120 times greater than that found in band 4. Therefore, when saturation occurs in bands 1, 2 and 3, the albedo of snow is estimated with the following equation:

In this equation, TM band 4 is used to represent the first part of the near-infrared region (0.725–1.00 μm). As in the visible portion of the spectrum, the selection of a representative band for the second portion of the near-infrared region (1.00–1.40 μm) presents an additional problem since no TM band falls in this region of the spectrum. It was found, through the inspection of several spectral reflectance curves, that snow reflectance in the 1.00–1.40 μm spectral region is on average 0.630 times less than that of TM band 4. Finally, TM band 7 is employed to represent reflectance in the mid-infrared (1.40–6.00 μm). Equations (4) and (5) are a modification ofReference Brest Brest's (1987) approach since TM bands are used instead of MSS. An improvement here is that a reflectance measurement is available for the mid-infrared portion of the spectrum (TM band 7).

In the model, the sensor is assumed to view the surface from directly overhead and the ground cover is assumed to be Lambertian with bidirectional reflectance distribution function (BRDF) f_r =ρ/π, where ρ is the surface reflectance. This is a reasonable assumption for freshly fallen snow in the visible wavelengths. However, as snow grain-size increases, the reflectance in the forward direction increases as well (Reference Hall, Chang, Foster, Benson and KovalickHall and others, 1989).

As a first approximation, surface reflectance of snow in each of the TM reflective bands is calculated as

where L[i] is the sensor radiance for a given TM band [i] associated with the surface, L _p[i] is the path radiance between the sensor and the surface at altitude z, T _v[i,z] is the vertical transmission from altitude z to the sensor, and E[i, z] is the total solar irradiance on an inclined surface at altitude z. In the calculation of reflectance, all components are computed from the radiative transfer model LOWTRAN7 for unobstructed horizontal surfaces (Reference KneizysKneizys and others, 1988).

The total in-band [i] solar irradiance for a sloping surface at altitude z is computed using Equation (equ.2). In contrast to the calculation of the total integrated shortwave irradiance (0.28–6.00μm), the last portion of the equation (0.5αE_h[i, z](l – coss)), which represents the contribution of terrain-reflected radiation, is not considered since a measurement of albedo needs to be known a priori (it is this variable that we are trying to estimate from remotely sensed measurements).

Surface thermal exitance

The total spectral radiance reaching a sensor in the thermal infrared can be expressed as (Reference Schott and VolchokSchott and Volchok, 1985; Reference Bartolucci, Chang, Anuta and GravesBartolucci and others, 1988)

where L[λ, z] is the total spectral radiance reaching the sensor at altitude z, ε_s is the emissivity of the surface observed, L _t[λ] is the black-body radiance associated with an object on the ground with temperature t, Τ_v[λ, z] is the spectral vertical transmission to altitude z, r is the reflectance of the surface observed (r=1-ϵ_s), L _D[λ,z] is the downwelled spectral atmospheric radiance, and L _U[λ,z] is the upwelled spectral atmospheric radiance (i.e. path radiance).

Surface-brightness temperature measurements are derived from Landsat TM thermal band 6 using the equation (Schott and Volchok, 1985):

where K ₁ is the first constant (60.776 m W cm⁻² sr⁻¹ μm⁻¹), Κ ₂ is the second constant (1260.56 K), and L[λ] is the spectral radiance reaching the sensor (m W cm ⁻² sr⁻¹ μm⁻¹). Atmospheric effects (i.e. path radiance and atmospheric transmission) can be removed by inverting Planck's equation and using an appropriate radiative transfer model. However, properly calibrated Landsat TM thermal infrared data, collected under clear atmospheric conditions (i.e. high visibility and low water vapour content), can be used to obtain accurate temperature measurements without the need for atmospheric corrections (Reference Bartolucci, Chang, Anuta and GravesBartolucci and others, 1988). This is because under clear skies, the atmospherically attenuated target radiance appears to be compensated by the path radiance.

Under clear sky conditions, surface thermal exitance (3.0–100.0 μm) can thus be determined through the Stefan-Boltzman relationship

All variables have previously been defined (see Equations (3) and (7)).

Procedure

In order to test the model a comparison is carried out with calculated and measured upward and downward fluxes for a cloudless atmosphere. The test area is Niwot Ridge, a long-term ecological research site, situated 30 km west of Boulder, Colorado, U.S.A. The surface of Niwot Ridge is a patchwork of snowfields and alpine tundra vegetation communities. In this area, strong prevailing winds from the west redistribute the snow. West-facing slopes and ridge tops are generally free of snow, while east-facing slopes usually accumulate a deep snowpack (4 m water equivalent or more). The deep snowfields are a very reliable source of runoff throughout the summer so a thorough understanding of their radiation balance income is desirable (Olyphant, 1986). On Niwot Ridge, variations in snow depth and the timing of snowmelt create vegetation gradients across topographic depressions.

For this experiment, data consisted of Landsat-5 Thematic Mapper (TM) data acquired under clear atmospheric conditions on 29 June 1984 (1712 h GMT), a 30 m resolution digital elevation model (DEM), and radiosonde data. The solar elevation and azimuth angles were 60° and 116°, respectively, on that date. Landsat TM bands 2, 4, and 7 were used to calculate albedo values of snowfields, and thermal band 6 of TM utilized in the determination of surface thermal exitance, and as an input in the calculation of longwave irradiance (Equation (equ.3)). The surface emissivity was specified as a constant (ε_s = 0.98) for all snow-covered surfaces.

Given estimates of downwelling radiation fluxes for unobstructed horizontal surfaces, at heights defined by the levels of the radiosonde ascents, the DEM was used to calculate values of shortwave and longwave flux densities, as well as path radiance and atmospheric transmission for surfaces of varying elevation, slope and aspect. Finally, surface albedo and thermal exitance values were estimated for the time of Landsat-5 overflight. The experiment was then carried a step further in order to simulate the radiation balance over the snowfields for the day (0600 to 1800h solar time). In that simulation, the snowpack was considered to be isothermal at 273 Κ with ε_s = 0.98, and surface albedo was assumed to vary between 0.50 and 0.56 as a function of the solar zenith angle, with the minimum at solar noon.

Table 1

Comparison of observed components of the radiation balance with model estimates at the time of satellite data acquisition (1712hGMT). Estimated values calculated over a 3 × 3 window over the snowfield in the Saddle area. Radiant flux densities are expressed in W m⁻²

Results

In order to evaluate the accuracy of the proposed method, the components of the radiation balance equation are compared with field measurements obtained during two summer field seasons (1984 and 1986) over a snowfield in the Saddle area of Niwot Ridge, and with data reported in the literature. Field observations consisted of surface albedo, surface brightness temperature, shortwave irradiance, and net all-wave radiation. The longwave irradiance was calculated as a remainder of the radiation-balance equation.

Estimated values of net radiation and its components for the snowfield are presented in Table 1. The predicted shortwave irradiance shows an overestimation of 7% compared with field measurements. Here, the differences between model predictions and observations should not only be interpreted as modelling errors. It is worth noting that the estimated value of the total shortwave irradiance is based on calculations over the spectral range 0.28–6.00 μm, a slightly broader spectral region than that covered by the Eppley pyranometer. This instrument measures the radiation within the 0.285–2.8 μm range. For longwave irradiance the model underestimates the observed value by a margin of 10%. Reference Lougeay and BrazelLougeay and Brazel (1982) have shown that models such as the one used to estimate the total longwave irradiance, tend to underestimate observed values. As stated by Lougeay and Brazel (1982, p. 234): “This may not be due to inefficiences of the model construct, but to advective or inversion effects on air temperature in the local areas.”

It can also be seen from Table 1 that the estimated albedo value is in error by 2% from the mean observed value. Some discrepancies between observed and estimated albedo value at this site may be related to differences in the state of the snowpack (contaminants, grain size, etc.) between dates at which data were recorded, and the isotropic assumption used by the model. The estimated value represents snow albedo in June, 1984 while the observed values correspond to field measurements in June, 1986. Albedo values of melting snow cover have been reported by Reference AguadoAguado (1985) in a study site in the Sierra Nevada for three consecutive field seasons. Aguado found mean albedo values of 0.49, 0.45 and 0.48 for the 1981, 1982 and 1983 melt seasons, respectively. Similarly, Reference LeDrew and WellerLeDrew and Weller (1978) measured albedo values of 0.48 to 0.57 over a decaying snow bank on Niwot Ridge. Both sets of values reported by LeDrew and Weller (1978) and Reference AguadoAguado (1985) are close to the value estimated from the model, and that observed during the 1986 field season on Niwot Ridge. The surface albedo map for the study area is presented in Figure Fig. 1.

An average brightness temperature value of 273 Κ was observed over the snowfield on Niwot Ridge (in 1986) using a PRT-10L infrared thermometer. A value of snow-surface temperature derived from the model over the test site was 272.63 K. This relatively small difference (0.37 K) is within TM band 6 radiometric sensitivity of 0.5 K. Translated into surface thermal exitance term, the model slightly underestimates the surface thermal exitance of the snowfield by 1.7 W m⁻². The map of surface thermal exitance derived with Landsat TM thermal band 6 can be found in Figure Fig. 2. In this figure, the details are not as sharp as in Figure Fig. 1 (albedo map) because data are recorded with a 120 m footprint in this portion of the spectrum.

Fig. 1.

Surface albedo map derived from Landsat-5 TM. Albedo values range from 0.06 to 0.60. Average albedo of snowfield in the Saddle area is 0.55.

Fig. 2.

Surface thermal exitance map derived from Landsat-5 TM thermal band 6. Thermal exitance values range from 276 to 500 W m⁻². Average thermal exitance of snowfield in the Saddle area is 309.68 W m⁻².

Finally, Figure Fig. 3 shows the results of the simulation of the components of the radiation budget over the snowfield on 29 June 1984 (0600–1800 h solar time). Field measurements for that date are not complete but field observations under similar clear sky conditions during that morning and other dates enable us to make some comparison. For example, a value of 906.68 W m⁻² was recorded for shortwave irradiance at 1000 h (solar time) at the Saddle weather station on Niwot Ridge. This station is situated only a few hundred metres from the snowfield. The curve of shortwave irradiance shows a value of 937.18 W m⁻², a difference of about 30 W m⁻². LeDrew and Weller (1978) have observed at the Saddle site a daily average (June–August) for longwave irradiance of 275 W m⁻², again a difference of approximately 30 W m⁻². The results presented herein are preliminary, and indeed indicate the need for further experimentation.

Fig. 3.

Simulated daily clear-sky shortwave irradiance (K↓), longwave irradiance (L↓), and net all-wave radiation (Q*) over a section of the snowfield in the Saddle area of Niwot Ridge (29 June 1984). The average slope of the surface is 4.76° and the aspect is 90.00°.