Anisotropic Reflection by Melting Glacier Ice

Anisotropic Reflection by Melting Glacier Ice: Measurements and Parametrizations in Landsat TM Bands 2 and 4 Wouter Greuell* and Martijn de Ruyter de Wildt*

T his article deals with the anisotropic reflection of radiation by melting glacier ice. Ground-based measurements of the directional distribution of the reflected radiation over the hemisphere (so-called BRDFs5bidirectional reflectance distribution functions) were made on the Morteratschgletscher (Switzerland) in Landsat TM bands 2 (520–600 nm) and 4 (760–900 nm). These BRDFs cover a wide range of solar zenith angles (26–758) and surface characteristics (quantified by a variation in the spectrally integrated albedo between 0.14 and 0.50). All BRDFs exhibit a similar pattern with a minimum in the nadir direction and a maximum in the forward limb, but the amount of anisotropy increases with increasing wavelength, with increasing solar zenith angle and with decreasing albedo. The data were used to derive parametrizations (one for each TM band) which relate the bidirectional reflectance (the reflectance in a specific direction) to the albedo for a given solar-view geometry. Specific parametrizations (one for each TM band) for “near-nadir reflection” are also presented. All these parametrizations can be used to convert satellite-derived bidirectional reflectances into surface albedos and thus to correct for anisotropic reflectance. The residual uncertainty in the albedo due to inaccuracy of the correction is estimated to be 0.02 in both TM bands. Elsevier Science Inc., 1999 INTRODUCTION In recent years, many studies have used satellite data to derive the surface albedo of glaciers (e.g., Hall et al., * Institute for Marine and Atmospheric Research (IMAU), Utrecht University, Utrecht, The Netherlands Address correspondence to Wouter Greuell, Inst. for Marine and Atmospheric Research (IMAU), Utrecht Univ., Princetonplein 5, NL 3584 CC Utrecht, The Netherlands. E-mail: [email protected] Received 30 January 1999; revised 11 May 1999. REMOTE SENS. ENVIRON. 70:265–277 (1999) Elsevier Science Inc., 1999 655 Avenue of the Americas, New York, NY 10010

1990; Winther, 1993; Knap and Oerlemans, 1996; Stroeve et al., 1997; Knap et al., 1999a; Reijmer et al., 1999). The motivation for these studies is a combination of two factors. The first is the fact that during the ablation season on glaciers absorption of short-wave radiation generally is the dominant source of energy. (see, e.g., Ambach, 1979; Greuell and Oerlemans, 1989; Duynkerke and van den Broeke, 1994; Hock and Holmgren, 1996). The absorbed fraction of the incoming short-wave radiation is, by definition, determined by the albedo, which may vary widely (Paterson, 1994). Therefore, the albedo exerts an enormous influence on the amount of melt. The second reason to study the surface albedo of glaciers by means of satellite data lies in the spatial variation of the albedo on many scales. Consequently, in many cases ground-based measurements are not representative (see, e.g., Knap and Oerlemans, 1996; Greuell et al., 1997; Stroeve et al., 1997). To summarize, the retrieval of surface albedos from satellite data is useful because ground-based measurements are often not representative for larger areas while the albedo exerts an enormous influence on the amount of melt. A problem of the “satellite-albedo studies” is the uncertainty introduced by the retrieval methods. Basically, there are four steps in the retrieval methods, all of which introduce uncertainty, namely, the calibration of the satellite sensors, the correction for atmospheric effects, the calculation of the spectrally integrated albedo from band albedos, and the correction for anisotropic reflection at the surface (see, e.g., Knap and Oerlemans, 1996). This article focuses on the correction for anisotropic reflection, which must be performed because satellite instruments measure radiances (radiance5rate of energy from a particular direction and contained within a unit solid angle; unit: W m22 sr21), whereas the surface albedo is the ratio of the (hemispheric) fluxes (also called “irradi0034-4257/99/$–see front matter PII S0034-4257(99)00043-7

266

Greuell and de Ruyter de Wildt

ances”) of the reflected and the incoming radiation. Therefore, we need to know how the outgoing flux is related to the radiance in the direction of the satellite or, in other words, how the albedo is related to the bidirectional reflectance (the reflectance in a specific direction) in the direction of the satellite. This can be done by measuring or modeling BRDFs (bidirectional reflectance distribution functions), which are functions describing how the reflectance of a surface varies with the solarview geometry. Stroeve et al. (1997) made a correction for anisotropic reflection based on modeling snow with a radiative transfer model. Model results agree well with field measurements over a small set of solar zenith angles. Corrections can rise to several tenths in percentage of the uncorrected albedo, depending on the solar-view geometry, wavelength, and snow grain size. The other “satellitealbedo studies” mentioned above used the assumption of isotropic reflection. However, Knap et al. (1999a) concluded that this is one of the least favorable assumptions used in the retrieval method. Similar conclusions were drawn by Winther (1993) and Reijmer et al. (1999). So corrections that should be made for the anisotropic reflection of solar radiation by glaciers are possibly significant. This article addresses the anisotropy of the reflected radiation field over melting glacier ice. Ground-based measurements relating to this topic were already made by Knap and Reijmer (1998). This resulted in five BRDFs of glacier ice obtained at five different sites with varying surface characteristics. However, each BRDF was collected for only one solar zenith angle, so the variation of the BRDF with the solar zenith angle is not well described by this data set. Because the solar zenith angle is expected to have a significant effect on the BRDFs, Knap and Reijmer (1998) were not able to deduce a general scheme for converting satellite-derived bidirectional reflectances into albedos. For the present study we performed ground-based BRDF measurements over a wide range of solar zenith angles and at three different locations with glacier ice of varying surface characteristics. Our aim was to cover as many as possible of the variations in BRDFs caused by variations in the solar zenith angle and the surface characteristics. As in the case of the measurements of Knap and Reijmer (1998), the measurements were carried out on the Morteratschgletscher (Switzerland) in Landsat TM Bands 2 (520–600 nm) and 4 (760–900 nm). These bands were chosen because knowledge of the albedos in these bands allows a fairly accurate estimate of the spectrally integrated albedo, as shown by Knap et al. (1999b). This so-called “narrow-to-broad band conversion” works due to the way the incoming flux and the albedo of glacier ice (see Cutler and Munro, 1996) vary with wavelength. Band 2 represents the relatively high albedos in the visible part of the spectrum, whereas Band 4 repre-

sents the transition to much lower albedos in the nearinfrared. At the same time these two bands, together with adjacent wavelengths, say 300–1000 nm, contain most of the incoming energy. Therefore, knowledge of the albedo for wavelengths larger than 1000 nm is not essential for an estimation of the spectrally integrated albedo. After introducing the coordinate system used to describe the solar-view geometry and some basic equations in the next section, we will describe the experimental set-up and present some measured BRDFs in the third section. The data will then be used to develop parametrizations that can be used to convert bidirectional reflectances into surface albedos (fourth section). Finally, in the fifth section, some attention will be paid to the importance that the derived corrections for anisotropic reflection have for the retrieval of the albedo of melting glacier ice. THE COORDINATE SYSTEM AND SOME BASIC EQUATIONS In this section the coordinate system (Fig. 1) and some basic equations are introduced. We follow the notation of Knap and Reijmer (1998). The surface albedo is defined as the ratio of the hemispheric fluxes of the reflected (Ir) and the incoming (Ii) radiation: I a5 r. Ii

(1)

The radiance (the intensity of radiation in a specific direction) R(hv, u) is related to Ir by Ir5#

0

0

# R(hv,u) sin hv cos hv dhv du. 2p p/2

(2)

The bidirectional reflectance r(hv, u) is defined as r(hv,u)5

R(hv,u) . Ii/p

(3)

For a Lambertian (5isotropically reflecting) surface the radiance and the bidirectional reflectance are independent of hv and u. In that case the bidirectional reflectance is equal to the albedo. Finally, the anisotropic reflection factor f(hv, φ) is defined as f(hv,u)5

r(hv,u) . a

(4)

Unless stated otherwise, all these variables refer to one of the TM bands. By combining Eqs. (1), (2), (3), and (4) it can be shown that p21#

0

#

0

2p p/2

f(hv,u) sin hv cos hv dhv du51,

(5)

or, in words, the anisotropic reflectance factor (f) has the

Anisotropic Reflection by Melting Glacier Ice

267

tive to each other after the field experiment by means of a Spectralon reflectance panel (see also Knap and Reijmer, 1998). An absolute calibration was not made. This would have been of no use for the present experiment because only the ratios of outgoing radiances to incoming fluxes were used for the analyses.

Figure 1. Coordinate system used to describe the solar-view geometry. The system consists of the solar zenith angle (hs), the view zenith angle (hv), and the relative azimuth angle (u). The latter is defined relative to the solar principal plane, that is, the sun is at u508. Consequently, forward scattering corresponds to u51808 and backscattering to u508. If the surface is occupied by surface roughness elements with a preferred orientation like sastrugi, the BRDF is also determined by the angle between the solar azimuth and the orientation of the roughness elements. In that case u has to be replaced by two variables, namely, the solar azimuth angle and the view azimuth angle (Warren et al., 1998). This was not done here, because the surfaces used for the measurements presented in this article did not show a preferred orientation of surface roughness elements. Note that all the angles are defined relative to the surface, which was never horizontal during the measurements discussed in this article.

property that its average value over the upward hemisphere, weighted by its contribution to the upward flux, is unity. THE MEASUREMENTS The Instruments The radiometers used for the present experiment were identical to those used by Knap and Reijmer (1998) and described in Knap (1997). The fluxes of the incoming radiation in the two Landsat TM bands were measured by two so-called Landsat global radiometers (LGRs). The reflected radiance in the same TM bands was measured with two Landsat direct radiometers (LDRs). The latter have a full opening angle of 5.060.28. In both bands the spectral response function of these radiometers is shifted by about 10 nm relative to the TM response functions (Knap and Reijmer, 1998), but since albedo and anisotropy vary slowly with wavelength, this wavelength shift causes only minor errors. Both LGRs and LDRs use silicon photo diodes with a response time of 1 ls to detect the radiation. The LDRs and LGRs were calibrated rela-

Procedure for Measuring Individual BRDFs In contrast to Knap and Reijmer (1998), who held the instruments by hand, we built a device to hold the LDRs and move them over the hemisphere (right-hand side of Fig. 2). This device resembles the “field-goniometer system” described in Sandmeier and Itten (1998). The surface viewed by the instruments is, for all instrument positions, situated in the center of the ring (diameter 134 cm). To obtain a BRDF, we always started to collect samples in the solar principal plane (variable hv, u508 and 1808). Sampling in a vertical plane was done by manually rotating the pseudovertical beam of the device around the point where it was connected to the ring on the surface, so that the sensors described a semicircle in the vertical plane. In fact, there are two rings on the surface. The inner ring is anchored to the ground; the outer ring, to which the upper part of the device is connected, can be moved over the inner ring. After completing the scan in the principal plane, the outer ring was moved (by hand) by 158 relative to the inner ring. Then, a second scan through a plane normal to the surface was made, etc., until after 12 scans the hemisphere was completely sampled. The measurement of a complete BRDF took about 5 min. During this time the LDRs sampled with a frequency of 0.5 Hz. The view zenith angle was also recorded at the same frequency by means of a potentiometer mounted at the point where the pseudovertical beam rotates around the ring. We estimate uncertainties in the angles describing the solar-view geometry (alignment errors) to be of the order of 2-38. The sampled view directions used to construct one of the BRDFs are given in Figure 3. Two kinds of samples are not plotted in the figure because they were discarded, namely, samples taken from the solar direction and samples for view zenith angles larger than 758, representing 7% of the reflected radiation for an isotropically reflecting surface. The latter were discarded because, at large view zenith angles, the sensors view the ring and the size of the viewed area increases dramatically. On average, individual BRDFs consist of 190 points. Note that the distribution of samples over the hemisphere is not representative for the distribution of the reflected radiation over the same hemisphere, for example, only 3% of the energy reflected by a Lambertian surface falls within 108 of the nadir direction. However, these directions are represented by, on average, 16% of our samples. So, directions within 108 from the nadir direction are oversampled. In order to compensate for this

268

Greuell and de Ruyter de Wildt

Figure 2. Experimental setup for the BRDF measurements. The incoming radiative fluxes in the two TM bands were measured with so-called Landsat global radiometers (LGRs). The third instrument adjacent to the LGRs is a pyranometer, an instrument measuring the incoming radiative flux over the entire solar spectrum. The radiances reflected by the surface in the two TM bands were measured with so-called Landsat direct radiometers (LDRs).

unrepresentative distribution of the measurements, the individual samples were weighted for the development of the parametrizations. This was done in such a way that the total weight of all the samples in each sector of 58 in the view zenith angle was proportional to the amount of radiation reflected by a Lambertian surface into all directions within that sector. The distance between the sensors and the surface in the center of the ring is roughly 125 cm. This means that the area viewed with the instruments in nadir position is a circle with a diameter of 11 cm. The area increases in size with the view zenith angle, for example, for the maximum view zenith angle of 758 it measures 11 cm (minor axis) by 43 cm (major axis). Hence, the surface viewed varies with the position of the sensor. It should be noted that this leads to errors if the surface is not homogeneous. The problem can be aggravated by misalignment of the experimental setup. During the entire sampling procedure described above, the incoming radiative fluxes in the two TM bands were measured with the LGRs at a frequency of 0.2 Hz (left-hand side of Fig. 2). The incoming fluxes were needed to calculate the bidirectional reflectance and the albedo of the surface (see the next section). This is the reason why the LGRs were not put in a horizontal position but placed parallel to the surface under investigation.

Finally, for each measurement location we determined the orientation of the surface in the middle of the ring. This information was needed to compute the solar zenith angle (defined with respect to the sloping surface). Figure 3. Sampled directions used for the construction of one of the BRDFs. The radial coordinate represents the view zenith angle, the azimuthal coordinate the relative azimuth direction in degrees (the backward direction is 08). The solar zenith angle during the measurement of this particular BRDF was 598.

Anisotropic Reflection by Melting Glacier Ice

Conditions during the Measurements As already stated in the introduction, all the measurements were carried out on the Morteratschgletscher (Switzerland) or, more precisely, on its main tributary, the Persgletscher (468249N, 98579E), at an elevation of approximately 2700 m a.s.l. During the experiment (6, 12, and 13 July 1998) the snow line on the Persgletscher was located at about 2900 m a.s.l., so the glacier surface was free of snow at the measurement locations. Retrieval of the surface albedo from satellite data is only possible in the absence of clouds. Therefore, one of the conditions stipulated for the measurements was that clouds should not obscure the sun. During the periods of the measurements no clouds or only a few clouds covered the sky, but, nevertheless, on some occasions a BRDF measurement had to be postponed because clouds obscured the sun. For the measurements three relatively plane and homogeneous surfaces without a preferred orientation of surface roughness elements were selected. The three surfaces had different surface characteristics. One surface was relatively white, the second surface relatively dark, and the third surface of intermediate brightness, as judged by eye. Although all three surfaces contained morainic material, we estimate that the concentrations were small enough for the BRDFs of these surfaces to be dominated by radiation reflected from the ice matrix and not by the morainic material. According to our subjective, rough estimate by eye 70–90% of the ice on the Persgletscher was rather similar to the ice of the investigated surface of intermediate brightness which, consequently, will be called “representative ice.” In terms of brightness, the other two surfaces represent rather extreme cases. During all of the measurements presented here, the three surfaces were melting. Their orientation was towards the north or northwest with slopes varying between 48 and 68. On each surface a series of measurements was carried out with, if cloud conditions were not adverse, a time interval of 30 min between individual BRDF measurements. This allowed us to investigate the effect that variation in the solar zenith angle had on the BRDFs (Fig. 4). The series on “white” and “representative ice” were started around solar noon and ended just before sunset, the series on “dark ice” started just after sunrise and ended around solar noon. For solar zenith angles larger than 508, that is, before 10:15 local time, the spectrally integrated (over the solar spectrum) albedo of the “dark ice” appeared to be similar in magnitude to the spectrally integrated albedo of the “representative ice” and even of the “white ice.” Around this time the albedo of the “dark ice” dropped quite dramatically. This is in agreement with our observation of the surface. During the first few hours after sunrise the investigated surface of the “dark ice” seemed

269

Figure 4. Conditions in terms of the solar zenith angle and the spectrally integrated albedo of the surface during the BRDF measurements presented in this article. During all the measurements the surface consisted of melting glacier ice. The spectrally integrated albedo was computed as a linear combination of the albedos in the two TM bands (asi50.427aTM210.354aTM4; see Knap et al., 1999b). The band albedos, in turn, were determined by integration over the hemisphere of the parametrizations of each measured BRDF [Eq. (8)].

to be rather similar to ice in the surrounding area, but later on it became much wetter and darker than the surrounding area. Therefore, the darkness and low albedo of the investigated surface at lower solar zenith angles are due not to high concentrations of morainic material but to accumulation of melt water. Most of the measured BRDFs were more or less symmetrical around the solar principal plane, an indication that generally the investigated surfaces were fairly homogeneous and that the alignment of the experimental setup was good enough. However, the first four BRDFs obtained over the “dark ice” showed considerably more asymmetry than the other measured BRDFs. These four BRDFs were excluded from the subsequent analysis, which will therefore be based on the remaining 38 BRDFs. The Resulting BRDFs To obtain individual BRDFs, measured outgoing radiances (R) were converted into bidirectional reflectances (r) by means of Eq. (3). Some characteristic BRDFs are plotted in Figure 5. The BRDF shown in panel a was obtained in TM band 2 over “representative ice” for an intermediate value of the solar zenith angle (hs5508). The three other panels show three effects on the BRDF, namely the effects of varying wavelength (panel b), solar zenith angle (panel c), and surface albedo (panel d). All the measured BRDFs, not only those shown in this figure, exhibit the same pattern. There is a minimum in the bidirectional reflectance for the nadir direction, or, in other words, limb brightening occurs for all azimuth directions. Limb

270

Greuell and de Ruyter de Wildt

Figure 5. BRDFs measured over melting glacier ice under various conditions: a) over ice of intermediate brightness, hs5508, in TM band 2; b) as a), but in TM band 4; c) as a), but hs5748; d) as a), but over “dark” ice. The band albedos were computed from the measured values of the bidirectional reflectance by integration of the parametrization of the individual BRDFs over the hemisphere [Eq. (8)]. In the plots the radial direction represents the view zenith angle (maximum 758) and the azimuth direction the relative azimuth angle. Contours give the bidirectional reflectance in increments of 0.04, except in the forward limb of panel c, where some contours are omitted.

brightening is most pronounced for the forward direction; it is weakest for the backward direction and therefore of intermediate strength for the sideward directions. The difference between the BRDFs is the amount of anisotropy as visualised in the figures by the contour density. The amount of anisotropy increases with increasing wavelength (panel b), with increasing solar zenith angle (panel c) and with decreasing surface albedo (panel d). The BRDF pattern revealed by our measurements is in broad agreement with the pattern found by Knap and Reijmer (1998). These authors also described the above-mentioned effects of wavelength and albedo. The effect of wavelength on the amount of anisotropy has got to do with the fact that the absorption coefficient of ice varies with wavelength (Grenfell and Perovich, 1981). Between 470 nm and 1030 nm the absorption coefficient increases monotonically with wavelength. So, in TM band 4 it is higher than in TM band 2. A higher absorption coefficient does not only lead to a lower albedo (compare Figure 5a with 5b) but also to a shorter mean path length through the ice of the photons emanating from the ice. A shorter path length means, on average, a smaller number of scattering events, which, in turn, leads to a better preservation of the angular distribution of the photons after a single scattering event, as described by the single scattering phase function. The single scattering phase function has a very strong peak

in the forward direction. In summary, the BRDF “remembers” more of the extreme asymmetry of the scattering between photons and individual grains if the wavelength is longer (Knap and Reijmer, 1988). Warren et al. (1998) explain why the amount of anisotropy increases with increasing solar zenith angle for snow. They state that “when the sun is high above the horizon, the forward scattered photons penetrate deeply into the snow and require many scattering events to be redirected upwards and eventually escape the snow pack. The emerging photons are therefore distributed more uniformly with angle into the upward hemisphere than is the case for lower solar elevations.” Variations in the surface albedo are caused by variations in the amount of morainic material and the amount of melt water on and near the surface. The effect of more morainic material and/or more melt water is similar to that of increasing wavelength: the absorption coefficient of the subsurface material increases, which leads to a lower albedo and a larger degree of anisotropy of the reflected radiation field. THE PARAMETRIZATIONS Introduction As stated in the introduction a major goal of this study is to develop parametrizations that describe the relation

Anisotropic Reflection by Melting Glacier Ice

between the bidirectional surface reflectance and the surface albedo for melting glacier ice. This will be done in this section for TM bands 2 and 4. From the measurements described in the previous section, it is obvious that the relation between the bidirectional reflectance and the albedo depends on the solar-view geometry, so the view zenith angle, the relative azimuth angle, and the solar zenith angle should appear as parameters in the equations. In the case where satellite data are processed, these angles can be computed for each pixel, provided that a topographic map with sufficient resolution is available. It is more of a problem to take account of the effect on the BRDFs of surface characteristics, like the amounts of morainic material and melt water at or near the surface. Satellite data do not give information about these surface characteristics. The only feasible solution to this problem is to describe the surface characteristics by the albedo itself. In order to make the development of general parametrizations, describing all the collected data in a single TM band, more transparent, we adopt a two-stage procedure. First, parametrizations of the individual BRDFs will be developed. This is done in order to study the part of the parametrizations relating to the view angles since solar zenith angle and surface characteristics (albedo) can be considered as constant during the measurement of a single BRDF (5 minutes). The coefficients in the parametrizations of each BRDF, as determined by regression analysis, are specific for each BRDFs. Therefore, they are a function of the solar zenith angle and the surface characteristics (albedo). In the second stage of the derivation of the general parametrizations, this functional relationship will be investigated. In principle, this leads in a straightforward way to the general parametrizations. However, in the case of melting glacier ice we will propose a simplification based on certain characteristics of the measured BRDFs. Parametrizations of Individual BRDFs Following Lindsay and Rothrock (1994) and Knap and Reijmer (1998), we anticipated that the parametrization of each BRDF would be of the following form: r5a01a1g11a2g21a3g3,

(6)

where g1, g2, and g3 are functions of the view angles hv and u, and are the same for all the BRDFs, whereas the coefficients ai differ for different BRDFs. Therefore, the ai vary with the solar zenith angle and the surface characteristics. For a certain set of functions gi the coefficients ai are determined with the criterion of minimum residual standard deviation in r. Knap and Reijmer (1998), following the example of Lindsay and Rothrock (1994), used the functions g15sin2 hv sin2 u, g25sin hv cos u and g35sin2 hv cos2 u. No arguments were given for their choice. We followed Lindsay and Rothrock (1994) and Knap and Reijmer (1998) by setting the number of functions gi equal to three, but we tried many different forms

271

of the functions, namely, gi5hi pi qi with hi and pi51, hv, sin hv, or cos hv and qi51, cos u or cos2 u. These forms of qi (and therefore of gi) imply symmetry around the solar principal plane. The best fit to the data, in terms of the residual standard deviation in r, was found for g15cos hv,

g25h2v cos u,

g35h2v cos2 u.

(7)

In that case the average residual standard deviation of the 38 analyzed BRDFs was 0.011 for TM Band 2 (on average, 88% of the variance explained) and 0.016 for TM Band 4 (on average, 89% of the variance explained). Using the above-mentioned functions of Lindsay and Rothrock (1994) and Knap and Reijmer (1998), the average residual standard deviation would have been 0.015 for TM Band 2 (85% of the variance explained) and 0.019 TM Band 4 (87% of the variance explained). The functions gi that provided the best fit to the data of the individual BRDFs [Eq. (7)] will be used in the next section to develop parametrizations that describe all the data in a single TM band in one equation. The parametrizations of the individual BRDFs were also used to compute the surface albedo by integrating r given by the parametrizations over the hemisphere [see Eqs. (4) and (5)]. a5p21#

0

#

0

2p p/2

r(hv,u) sin hv cos hv dhv du.

(8)

One Parameterization for All BRDFs Given the parametrizations of the individual BRDFs [Eqs. (6) and (7)], it is straightforward to develop a parametrization that describes all the data of a single band in one equation. The coefficients ai in Eq. (6) are independent of the view angles, but they differ for different BRDFs. Hence, variations in ai are due to variations in the solar zenith angle (hs) and the surface characteristics (a). Introducing the coefficients bi5ai/a yields for each BRDF [Eq. (9)]: f5b01b1g11b2g21b3g3,

(9)

where the bi again depend on the solar zenith angle and the albedo. Inserting this equation for f into the normalization condition given by Eq. (5) and with the functions gi proposed in the previous section [Eq. (7)], we obtain

1

2

1

2

2 1 1 f511b1 g12 1b2g21b3 g32 p21 . 3 16 4

(10)

It is important to notice that the values of the constants appearing in this equation depend on the choice of the functions gi. It should also be noted that, irrespective of how the coefficients bi vary with the solar zenith angle and the albedo, this equation ensures the validity of Eq. (5). To find the functions bi(hs, a) for melting glacier ice, it would have been straightforward to insert the values of ai and a found for the individual BRDFs into bi5ai/

272

Greuell and de Ruyter de Wildt

Figure 6. Amount of anisotropy, expressed as the standard deviation in the anisotropic reflectance factor, against solar zenith angle for all the BRDFs obtained over “white ice”. The albedo remained within a narrow range (0.62-0.67 in TM band 2 and 0.53-0.60 in TM band 4). The two lines show exponential fits to the points.

a, giving 38 values of b1, b2, b3, hs, and a for each TM band. This data set could then be used to find optimal functions bi(hs, a). In fact, we followed this procedure, but the result, in terms of the residual standard deviation in f, was no better than the result of an alternative procedure. This alternative leads to a more transparent equation with fewer coefficients. It uses a specific property of the measured glacier-ice BRDFs, namely, the fact that the pattern (minimum for the nadir direction, etc.) remains almost the same for varying solar zenith angles and albedos. If the pattern is independent of hs and a, Equation (10) can be simplified as follows:

31

2

1

2 1 f511 c1 g12 1c2g21c3 g32 p2 3 16

24

1 1 Z(hs)A(a) 4

(11)

where the ci are constants that are independent of hv, u, hs and a, and Z and A are functions of hs and a, respectively. This form of the parametrization ensures that effects of variations in hs and a on the three functions of hv and u between square brackets are identical, which means that the pattern does not vary with hs and a. From the measurements (Fig. 5) it appeared that the amount of anisotropy varies with hs and a. This effect is described by the reinforcement factors Z(hs) and A(a). We determined the form of these factors by first calculating the standard deviation of f for each BRDF. Then, this quantity was plotted against the solar zenith angle for a small range of the albedo (Fig. 6) and against the albedo for a small range of the solar zenith angle (Fig. 7). The assumption that the reinforcement factors are proportional to the standard deviation in f leads to Z(hs)5exp(hs/hc),

(12)

A(a)5a2m.

(13)

In order to obtain an optimal fit to all the data the coef-

ficients hc and m were still considered as free parameters. To summarize, it is suggested that glacier-ice BRDFs obey Eqs. (11), (7), (12), and (13). The unknown constants are c1, c2, c3, hc and m. These were found by performing a regression analysis on all the individual samples using the criterion of minimum residual standard deviation in f. The results are given in Table 1. In both TM bands, the parametrizations explain 89% of the variance. The largest contribution to the explained variance comes from g1, the smallest contribution from g3. The total explained variance is reduced by 3% (TM band 2) and 2% (TM band 4) if m is set equal to 1. This result is also shown in the table, since the calculation of a from r is easier for m51 than for m?1 (see subsection after next). For one situation (TM band 2, hs5508, band al-

Figure 7. Amount of anisotropy, expressed as the standard deviation in the anisotropic reflectance factor, against band albedo for all the BRDFs obtained while the solar zenith angle ranged between 268 and 368. The two lines show power fits to the points of each of the two TM bands.

Anisotropic Reflection by Melting Glacier Ice

273

Table 1. Results of the Regression Analyses Used To Fit Eq. (11) to the Dataa TM Band

Residual St. Dev.

Variance Explained

By g1

By g2

By g3

c1

c2 (rad22)

c3 (rad22)

hc (degs)

m

2 2 4 4

0.036 0.043 0.048 0.053

89% 86% 89% 87%

56% 53% 58% 56%

21% 22% 22% 22%

13% 11% 9% 9%

20.0123 20.0292 20.0284 20.0509

20.00374 20.00810 20.00726 20.01280

0.00225 0.00462 0.00314 0.00562

25 30 30 36

1.7 1.0 1.4 1.0

a Given are the residual standard deviation in f, the amount of variance explained by the entire parametrization, the contributions made to the explained variance by the individual functions gi and the optimal values of the coefficients. For each TM band two results are provided. The first was obtained with m as a free parameter; the second was obtained for m51.

bedo50.52) the BRDF calculated by the parametrization is compared to the measurements in Figure 8. For this particular BRDF the residual standard deviation is 0.013. Discussion of the discrepancies is postponed until the subsection after next. Parametrization of “Near-Nadir Reflectance” The parametrizations proposed in the previous section provide the best fit to all the measured anisotropic reflectance factors which are distributed over almost the entire hemisphere (Fig. 3). However, over terrain with small slopes, Landsat TM sensors “look” only into “nearnadir directions.” Over horizontal terrain the view zenith angle of these sensors has a maximum of 7.78 (Lillesand and Kiefer, 1987, p. 568). For these “near-nadir view angles” two specific parametrizations (for the two TM bands) were developed by means of the measurements presented in this article. There was a twofold motivation to do this. First, for “near-nadir directions” specific “near-nadir parametrizations” will be more accurate than parametrizations valid for the entire hemisphere. Second, the parametrizations can be simplified. This is due to the fact that variations in f for “near-nadir directions” are small (see Fig. 5). Owing to these small variations, the part of the parametrization describing the view-angle dependence [the part between square brackets in Eq. (11)] can be skipped for “near nadir” view angles.

Figure 8. Comparison of one of the measured BRDFs with the BRDF computed by means of the general parametrization [Eq. (11)] for the same conditions, that is, TM band 2, hs5508, band albedo50.52. The measured BRDF was also plotted in the upper left panel of Figure 5.

For the development of the “near-nadir parametrization” all the measurements obtained within 108 of the normal to the surface were selected. This yielded, on average, 31 “near-nadir samples” per BRDF. For each BRDF the mean f of these samples was computed and called the “near-nadir anisotropic reflectance factor” (fn). The fn were then plotted against hs and a (as in Figs. 6 and 7) in order to determine the form of the relationships between fn, on the one hand, and hs and a, on the other hand. Whereas the entire BRDF was found to be represented optimally by an exponential relationship between f and hs [Eq. (12)], the “near-nadir data” are best represented by a linear relationship between fn and hs. The relation between fn and a is again best described by a power law, so fn can be written as fn511(c41c5hS)a2mn.

(14)

For both TM bands the values of the constants c4, c5, and mn were determined by fitting Eq. (14) to the data sets of fn, hs, and a. Results are given in Table 2. The explained variances in fn are extremely high (99% for TM band 2 and 98% for TM band 4). Inversion of the Problem If the solar-view geometry and the surface albedo are known, the bidirectional reflectance can be calculated directly with the parametrizations given by Eqs. (11) or (14). However, in satellite retrieval one faces the inverse

274

Greuell and de Ruyter de Wildt

Table 2. Results of the Regression Analyses Used To Fit Eq. (14) to the “Near-Nadir Data”a TM Band

Residual St. Dev.

Variance Explained

c4

c5 (rad21)

mn

2 4

0.008 0.013

99% 98%

0.02553 0.01496

20.08708 20.09786

1.6 1.3

a Given are the residual standard deviation in fn, the amount of variance explained by the parameterization and the optimal values of the coefficients.

problem: The albedo must be computed from a value of the bidirectional reflectance. To solve this problem Eq. (11) [or Eq. (14)] is rewritten as follows: r f5 511f˜a2m, a

(15)

where f and m can be replaced by fn and mn, respectively, in the case of “near-nadir calculations.” The symbol f˜ represents the factors in Eq. (11) [or Eq. (14)] which describe the effect of the solar-view geometry. For the calculation of the albedo from the bidirectional reflectance Eq. (15) is rearranged as follows: a5r2f˜a12m (16) This equation can only be solved analytically for some integer values of m. If m51, the difference between a and r (“the correction to r”) is independent of the albedo, as suggested by Knap and Reijmer (1998). However, optimum values for m (and mn) are between 1 and 2. In that case, Eq. (16) must be solved numerically. This was done for all our measurements of r (7227 samples) by means of a Newton algorithm (Atkinson, 1978), where f˜ was calculated with Eq. (11). Most of the solutions (99.75%) were found within five iterations. The residual standard deviations in the computed albedos were found to be 0.018 in TM band 2 and 0.020 in TM band 4. The same values of the residual standard deviations were found for m51. Assuming isotropic reflectance, the residual standard deviations would have been 0.058 (TM band 2) and 0.067 (TM band 4). The success of the parametrizations in estimating the albedo from values of the bidirectional reflectance depends on the solar-view geometry. Figure 9 shows the residual standard deviation as a function of the view geometry. For most angles the residual standard deviation is less than 0.02, but it increases dramatically with view zenith angle in the forward limb. Figure 10 shows how the residual standard deviation varies with the solar zenith angle. Each symbol represents the residual standard deviation for one of the BRDFs. For solar zenith angles less than 668 and for “representative” and “white ice” the residual standard deviation remains less than 0.015. However, for the same kinds of ice, the success of the parametrization diminishes rapidly for larger solar zenith angles. The parametrization is also less accurate for the “dark ice.” In summary, it appears that corrections for

anisotropic reflection by means of the proposed parametrizations have an accuracy better than approximately 0.02, except for large view zenith angles in the forward limb and for large solar zenith angles (.668). The accuracy also decreases with decreasing albedo. CONCLUSIONS AND DISCUSSION In this article BRDF measurements of melting glacier ice in two Landsat TM bands are presented. The measurements were performed on the Morteratschgletscher in Switzerland. In total, for both TM bands 2 and 4, 38 BRDFs were collected, over three surfaces with varying characteristics (quantified by a variation in the spectrally integrated albedos between 0.14 and 0.50) and under varying solar zenith angles (ranging between 268 and 758). The temporal variation in the albedo over the third, “dark,” surface indicates that a large part of these albedo

Figure 9. Residual standard deviation in the computed band albedo (TM2) as a function of the view angles (coordinates as in Fig. 5). Calculated albedos were obtained from all the measured bidirectional reflectances by means of the general parametrization [Eq. (11)]. Measured albedos were obtained by integration of the parametrizations of the individual BRDFs [Eq. (8)]. Subsequently, the residual standard deviation was computed in bins of 7.58 for the view zenith angle and 158 for the relative azimuth angle.

Anisotropic Reflection by Melting Glacier Ice

Figure 10. Residual standard deviation in the computed band albedo (TM2) as a function of the solar zenith angle for three types of melting glacier ice. Each symbol represents a single BRDF.

differences was caused by variations in the amount of water at or near the surface. All BRDFs exhibit the same pattern. There is a minimum in the bidirectional reflectance for the nadir direction. Therefore, limb brightening occurs for all azimuth directions, but it is most pronounced for the forward direction and is weakest for the backward direction. Although the pattern is almost invariable, the amount of anisotropy increases with increasing wavelength, with increasing solar zenith angle and with decreasing surface albedo. For each TM band two parametrizations relating the bidirectional reflectance to the albedo were developed. The first describes the bidirectional reflectance for all directions of the upper hemisphere [Eq. (11)], with the functions gi given by Eq. (7), the functions Z and A given by Eqs. (12) and (13), respectively, and the coefficients ci, hc, and m given in Table 1. The second parametrization [Eq. (14)] was developed with the measurements from “near-nadir angles” (view zenith angle less than 108) only. Therefore, it is more accurate for these angles but not valid for larger view zenith angles. This means that it is only applicable to Landsat TM data for parts of glaciers or ice sheets with relatively gentle slopes. The values of the coefficients c4, c5, and mn in Eq. (14) are given in Table 1. We suggested that the inverse problem of computing the albedo from the bidirectional reflectance could be solved with a Newton algorithm. For “near-nadir viewing directions” and for TM band 2, the correction for anisotropic reflectance, which is the difference between the albedo and the bidirectional reflectance, is shown in Figure 11 as a function of the solar zenith angle and the bidirectional reflectance itself. Remember that the albedo is set equal to the bidirectional reflectance under the assumption of isotropic reflection. The correction increases with increasing solar

275

Figure 11. Correction for anisotropic reflection for “near-nadir directions” in TM band 2 as a function of the solar zenith angle and the bidirectional reflectance [Eq. (14)]. In order to obtain the albedo, the correction must be added to the bidirectional reflectance. Values derived from two Landsat images of Breidamerkurjo¨kull (outlet of Vatnajo¨kull, Iceland) correspond to the double arrow. The ranges in bidirectional reflectance and solar zenith angle shown correspond to those of our measurements.

zenith angle and decreasing bidirectional reflectance because the amount of anisotropy varies in the same way. To give an example, bidirectional reflectances in TM band 2 of melting glacier ice on Breidamerkurjo¨kull (outlet of Vatnajo¨kull, Iceland) ranged between 0.16 and 0.27, as derived from two Landsat images (see also Reijmer et al., 1999). At the time these images were collected, the solar zenith angle was 548. Therefore, the correction for anisotropic reflection in TM band 2 ranges between 0.10 and 0.12. The corresponding corrections in TM band 4 range between 0.11 and 0.12, yielding a spectrally integrated correction (equation in caption of Fig. 4) between 0.08 and 0.09. Obviously, this is a significant correction, so the assumption of isotropic reflection may lead to large errors in satellite-derived albedos. Indeed, application of this correction to the two abovementioned Landsat images leads to significantly better agreement between ground-based measurements of the albedo and satellite-derived albedos (Carleen Reijmer, personal communication). If the parametrizations for the entire hemisphere [Eq. (11)] are used to calculate the albedos from the measured bidirectional reflectances, the residual standard deviation in the albedo is 0.018 in TM band 2 and 0.020 in TM band 4. These numbers can be taken as a measure of the uncertainty in the band albedos due to the uncertainty in the correction for anisotropic reflection [Eq. (11)]. However, there are two reasons for ad-

276

Greuell and de Ruyter de Wildt

justing these uncertainty estimates. Neither can be quantified. Calculated residual errors are probably caused mainly by inhomogeneities of the observed surface and by alignment errors. These errors will be partly random. Therefore, the given residual standard deviation is an overestimation of the uncertainty. On the other hand, the measurements were obtained at only three locations. Although the choice of these locations was based on the representativity of the investigated surfaces for all the ice on the glacier, there will certainly be glacier-ice surfaces which are badly described by the parametrizations derived from the data obtained at the three selected locations. Therefore, the given residual standard deviation is an underestimation of the uncertainty. In principle, the proposed parametrizations can be used without any other knowledge than that of the solarview geometry and the assurance that the surface consists of melting glacier ice. The parametrizations should not be applied to areas where the amount of morainic material is so abundant that a significant part of the reflected radiation emanates from the morainic material, nor should they be applied to snow. The reason is that the BRDFs of morainic material and of snow can differ significantly from those of glacier ice. For snow this appears from BRDF measurements that we performed with the experimental setup described in this article. Although these BRDFs are of lesser quality than the ice BRDFs, they are good enough to reveal a completely different pattern. For smaller solar zenith angles the minimum in the reflectance lies in the limb of the backward radiation. The minimum shifts towards nadir for larger solar zenith angles. A similar picture emerges from snow BRDFs obtained at the top of the atmosphere (Suttles et al., 1988). The parametrizations do not apply outside the ranges in band albedo and solar zenith angle prevailing during our measurements. Extrapolation beyond these ranges hinges on the correctness of the form of the factors describing the effect of variations in albedo and solar zenith angle [Eqs. (12), (13), and (15)]. Possibly simulations of BRDFs by means of a radiative transfer model of ice [see, e.g., Choudhury and Chang (1981) for snow] could determine whether these equations have a physical basis or not. This is also true for the three functions that describe the variations in reflectance with the view angles [Eq. (7)]. Finally, we would like to discuss whether the presented parametrizations for the two TM bands could perhaps be applied to data collected by different satellite sensors, representing different wavelengths. This discussion is guided by the almost coincidence in Figure 7 of the two curves representing the two TM bands. To understand the implications of this “almost coincidence,” two hypothetical BRDFs are compared, characterized by the same band albedo. BRDF 1 is obtained in TM band 2, BRDF 2 in TM band 4. Because wavelength affects

the band albedo, the band albedos of the two BRDFs can only be the same if the surficial layer of BRDF 1 contains more absorbing material than the surficial layer of BRDF 2. The “almost coincidence” of the two curves in Figure 7 means that BRDF 1 and 2 have almost the same amount of anisotropy and are therefore almost identical, despite the differences in wavelength and the amount of absorbing material. We conclude that the band albedo appearing in Eq. (13) and on the right-hand side of Eq. (14) represents the effects on the BRDFs of both the absorbing material and of wavelength. Qualitatively this can be understood because absorbing material in and on the ice and wavelength both affect the radiation field by changing the absorption coefficient (see the third section). In fact, the two parametrizations for the two TM bands are almost identical. The practical implication is that the presented parametrizations can be applied at other wavelengths near to those of TM bands 2 and 4. The experimental setup conceived and constructed by Henk Snellen, Peter Braam, and Wim Boot was a success, although it was used for the first time. Wouter Knap, Alexander Los, Carleen Reijmer, and Theo Gerkema helped us with various parts of the work that led to this article. An earlier version of this article was reviewed by Wouter Knap, Richard Bintanja, and Roderik van de Wal. Sheila McNab corrected the English. We are very grateful to all of the above-mentioned persons. This study was financed by NOP-II Project 013 001236.10, which is a part of the Dutch National Research Programme on Global Air Pollution and Climate Change, and by Project EO030, “Glacier Mass Balance from Space,” financed by SRON (Space Research Organisation Netherlands).

REFERENCES Ambach, W. (1979), Zum Wa¨rmehaushalt des Gro¨nla¨ndischen Inlandeises: Vergleichende Studie im Akkumulations- und Ablationsgebiet. Polarforschung 49(1):44–54. Atkinson, K. E. (1978), An Introduction to Numerical Analysis, Wiley, New York, 587 pp. Choudhury, B. J., and Chang, A. T. C. (1981), On the angular variation of solar reflectance of snow. J. Geophys. Res. 86(C1):465–472. Cutler, P. M., and Munro, D. S. (1996), Visible and near-infrared reflectivity during the ablation period on Peyto Glacier, Alberta, Canada. J. Glaciol. 42:333–340. Duynkerke, P. G., and van den Broeke, M. R. (1994), Surface energy balance and katabatic flow over glacier and tundra during GIMEx-91. Global Planet. Change 9:17–28. Grenfell, T. C., and Perovich, D. K. (1981), Radiation absorption coefficients of polycrystalline ice from 400–1400 nm. J. Geophys. Res. 86 (C8):7447–7450. Greuell, W., and Oerlemans, J. (1989), Energy balance calculations on and near Hintereisferner (Austria) and an estimate of the effect of greenhouse warming on ablation. In Proceedings of the Symposium on Glacier Fluctuations and Climatic Change in Amsterdam (J. Oerlemans, Ed.), Kluwer Academic, Dordrecht, pp. 305–323.

Anisotropic Reflection by Melting Glacier Ice

Greuell, W., Knap, W. H., and Smeets, P. C. (1997), Elevational changes in meteorological variables along a midlatitude glacier during summer. J. Geophys. Res. 102(D22): 25,941–25,954. Hall, D. K., Bindschadler, R. A., Foster, J. L., Chang, A. T. C., and Siddalingaiah, H. (1990), Comparison of in situ and satellite-derived reflectances of Forbindels Glacier, Greenland. Int. J. Remote Sens. 11:493–504. Hock, R., and Holmgren, B. (1996), Some aspects of energy balance and ablation of Storglacia¨ren, northern Sweden. Geogr. Ann. 78A (2–3):121–132. Knap, W. H. (1997), Satellite-derived and ground-based measurements of the surface albedo of glaciers. Ph.D. thesis, Institute of Marine and Atmospheric Sciences, Rijksuniversiteit Utrecht, 175 pp. Knap, W. H., and Oerlemans, J. (1996), The surface albedo of the Greenland ice sheet: satellite-derived and in situ measurements in the Søndre Strømfjord area during the 1991 melt season. J. Glaciol. 42:364–374. Knap, W. H., and Reijmer, C. H. (1998), Anisotropy of the reflected radiation field over melting glacier ice: measurements in Landsat TM bands 2 and 4. Remote Sens. Environ. 65:93–104. Knap, W. H., Brock, B. W., Oerlemans, J., and Willis, I. C. (1999a), Comparison of Landsat-TM derived and groundbased albedos of Haut Glacier d’Arolla, Switzerland. Int. J. Remote Sens., in press.

277

Knap, W. H., Reijmer, C. H., and Oerlemans, J. (1999b), Narrowband to broadband conversion of Landsat TM glacier albedos. Int. J. Remote Sens. 20 (10):2091–2110. Lillesand, T. M., and Kiefer, R. W. (1987), Remote Sensing and Image Interpretation, Wiley, New York, 721 pp. Lindsay, R. W., and Rothrock, D. A. (1994), Arctic sea ice albedo from AVHRR. J. Clim. 7:1737–1749. Paterson, W. S. B. (1994), The Physics of Glaciers, Pergamon, Oxford, 380 pp. Reijmer, C. H., Knap, W. H., and Oerlemans, J. (1999), The surface albedo of the Vatnajo¨kull ice cap, Iceland: a comparison between satellite-derived and ground-based measurements. Boundary-Layer Meteorol. 92 (1):125–144. Sandmeier, St., and Itten, K. I. (1999), A field goniometer system (FIGOS) for acquisition of hyperspectral BRDF data. IEEE Trans Geosci. Remote Sens. 37 (2):978–986. Stroeve, J., Nolin, A., and Steffen, K. (1997), Comparison of AVHRR-derived and in situ surface albedo over the Greenland ice sheet. Remote Sens. Environ. 62:262–276. Suttles, J. T., Green, R .N., Minnis, P., et al. (1988), Angular Radiation Models for Earth–Atmosphere System. Volume I: Shortwave Radiation, NASA Reference Publication 1184, NASA, Washington, DC. Warren, S.G., Brandt, R. E., and O’Rawe Hinton, P. (1998), Effect of surface roughness on bidirectional reflectance of Antarctic snow. J. Geophys. Res. 103(E11):25,789–25,807. Winther, J.-G. (1993), Landsat TM derived and in situ summer reflectance of glaciers in Svalbard. Polar Res. 12:37–55.