Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography

. Accurate knowledge of temperate glacier mass balance is essential to understand the relationship between glacier and climate. Deﬁned as the reﬂected fraction of incident radiation over the whole solar spectrum, the surface broadband albedo is one of the most important variable in a glacier’s mass balance. This study presents a new method to retrieve the albedo of frozen surfaces from terrestrial photography at visible and near infrared wavelengths. This method accounts for the anisotropic reﬂectance of snow and ice surfaces and uses a radiative transfer model for narrow-to-broadband conversion. The accuracy of the method was assessed using concomitant measurements of albedo during the summers 2008 and 2009 on Saint Sorlin Glacier (Grandes Rousses, France). These albedo measurements are performed at two locations on the glacier, one in the ablation area and the other in the accumulation zone, with a net radiometer Kipp and Zonen CNR1. The main sources of uncertainty are associated with the presence of high clouds and the georeferencing


Introduction
The energy transfer at the surface of a glacier is controlled by meteorological variables (e.g., air temperature and humidity profiles, wind speed) as well as surface variables (e.g., albedo, surface temperature, surface roughness) (Hock, 2005). Temperate glaciers, which are characterized by a temperature being close to the melting point (Paterson, 1994, p. 12), are very sensitive to these variables (Vincent, 2002).
Correspondence to: M. Dumont (marie.dumont@meteo.fr) The broadband bihemispherical reflectance of a surface (Bonnefoy, 2001) -referred to as albedo in the following -is defined as the fraction of the incident irradiance that is reflected by the surface over the whole solar spectrum. Thus, the albedo governs the amount of shortwave radiation that is effectively absorbed by the material. Since shortwave radiation plays a major role in the energy budget of temperate glaciers (Sicart et al., 2008;Slaymaker and Kelly, 2007, p. 57), the albedo is one of the leading variables controlling the energy balance . Consequently, an accurate estimation of the surface albedo is essential to compute correctly the mass balance of glaciers.
The albedo of glaciers is highly variable, both temporally and spatially. It can range from more than 90 % for fresh snow down to around 20 % for dirty ice (Warren, 1982;Oerlemans and Knap, 1998). These variations are due principally to the variability in size and shape of the ice grains, the distribution, type, and content of impurities, the surface roughness, the liquid water content, and the characteristics of the irradiance (Warren, 1982;Aoki et al., 2007;Jin et al., 2008;Dozier et al., 2009). In the context of modeling the mass balance of glaciers and/or modeling global climate, spatial and temporal estimates of the albedo of ice and snow are required. In practice, this variable is often parameterized (Hock, 2005;Gardner and Sharp, 2010). However, albedo values obtained with such empirical methods often exhibit substantial discrepancies with ground measurements (Pedersen and Winther, 2005). Therefore, more systematic measurements are needed to develop, validate, and improve parameterizations schemes.
Terrestrial and spaceborne remote sensing has proven to be a desirable and efficient way to monitor the temporal and spatial variations of snow properties (Corripio, 2004;Dozier et al., 2009). For instance, studies have reported on the use of Published by Copernicus Publications on behalf of the European Geosciences Union. remote sensing to retrieve snow albedo (Dozier et al., 2009), fraction of snow covered area , or grain size (Lyapustin et al., 2009;Jin et al., 2008). The retrieval of snow and ice albedo from remote sensing techniques is based on the combination of measured spectral reflectance which often implies errors larger than 10 % (Dozier et al., 2009).
In addition it has been shown that ignoring the anisotropic reflection pattern of snow surfaces in the computation of albedo can introduce an additional error on the order of 10 %, especially when the incident solar zenith angle is high (Painter and Dozier, 2004;Hudson et al., 2006;Dumont et al., 2010). Most of these studies propose different solutions to address the two main difficulties encountered while converting remote sensing reflectance measurements into albedo values: (i) the fact that snow and ice are not Lambertian surfaces, meaning that the reflected radiance varies over the overlying hemisphere (Hudson et al., 2006(Hudson et al., , 2010Lyapustin et al., 2009); (ii) the conversion from narrow band measurements to broadband values required when computing the radiative budget of the surface. Corripio (2004) developed a method to retrieve glacier albedo from terrestrial photography that addressed the topographic effects on incident and reflected radiation. Nevertheless, this method does not take into account the anisotropy of the surface. In addition, it requires at least one measurement of albedo to be made in the field at a location visible in the image to estimate the albedo in the remaining pixels. This shortcoming arises from the fact that the system designed by Corripio (2004) is not readily capable of providing estimates of the spectral ground reflectance to allow a direct spectral conversion of the image into broadband albedo. In addressing most topographic and atmospheric effects on the irradiance, Sirguey et al. (2009) proposed a method to retrieve spectral ground reflectance values in mountainous terrain from remotely sensed data. Although this technique did not readily implement the effect of directional reflectance, it provided an incremental step towards the direct mapping of the surface albedo in rugged terrain.
Building on the work of Corripio (2004) and Sirguey et al. (2009), this study reports on the design of a new method to retrieve glacier surface albedo from terrestrial photography, while accounting for the anisotropic reflectance of the frozen surface. In the first section, a precise definition of the albedo and the related radiometric quantities used in the paper are provided. The study site (Saint Sorlin Glacier, French Alps) and measurement devices are presented in Sect. 3. Section 4 describes the algorithm that was designed to map the albedo from visible and near IR digital photography of the glacier. The results obtained with this method on the validation site during the 2008 and 2009 summers are presented in Sect. 5. Finally, the advantages and accuracy of the method are discussed in Sect. 6.

Definitions
The Bidirectional Reflectance Distribution Function (BRDF) -Strub et al., 2006;Nicodemus et al., 1977). Thus, the BRDF is formulated as where λ is the wavelength. θ i and θ v (∈ [0,90 • ]) are the zenith angles of the incident and viewing directions, respectively, φ i and φ v (∈ [0,360 • ]) are the corresponding azimuth angles. The relative geometry of the source, the target, and the sensor is illustrated in Fig. 1. In general, the relative azimuth |φ v −φ i | is sufficient to describe the BRDF (Warren et al., 1998). In addition, when there is no preferred orientation of the surface roughness in the azimuth direction, the BRDF can be considered to be symmetric with respect to the solar principal plane. The latter contains the incident beam and the normal to the surface (Hudson et al., 2006). Accordingly, the BRDF can be assumed to depend only on the absolute value of the relative azimuth, The spectral albedo α(θ i ,λ), also referred to as the directional-hemispherical reflectance (Nicodemus et al., 1977;Schaepman-Strub et al., 2006), is the ratio between the reflected radiant flux integrated over the upper hemisphere and the incident collimated flux as a function of the wavelength. Using Eq. (1), the spectral albedo can be written In order to isolate the angular dependence of the snow surface reflection, the anisotropic reflectance factor R(θ i ,θ v ,φ,λ) can be formed (Hudson et al., 2006). It measures the fraction of the radiant flux effectively reflected by the surface to that which would be reflected if the surface was perfectly diffuse (i.e., lambertian) under the same condition of illumination. Thus, R is simply a normalization of the BRDF value given by In the case of lambertian surfaces R(θ i ,θ v ,φ,λ) is uniformly equal to 1. The albedo represents the reflected fraction of the total hemispherical incident flux integrated over the whole solar spectrum λ (Bonnefoy, 2001). It identifies to the bihemispherical reflectance defined by Schaepman-Strub et al. (2006) and integrated over λ. It can be obtained by integration of the spectral albedo given in Eq. (2) weighted by the incident spectral solar radiance L i (θ i ,φ i ,λ) originating  from an elementary solid angle in direction (θ i ,φ i ). Thus, the albedo is given by In the following, the subscript v refers to the viewing direction of the camera and s to the sun. Tilded angles are defined in the reference frame attached to the considered pixel of the glacier (see Fig. 1b and Appendix A).

Saint Sorlin Glacier
Situated by 45.10 • N and 6.10 • E in the Western Alps of France (Grandes Rousses area), Saint Sorlin Glacier covers 3 km 2 (Fig. 2). The glacier extends from around 2700 m a.s.l. at its terminus to nearly 3500 m a.s.l. at Etendard peak. Its mass balance has been monitored by the Laboratoire de Glaciologie et de Geophysique de l'Environnement (LGGE, Grenoble, France) for several decades (the monitoring program started in 1957) (Vincent, 2002). Much information about this site is available at http://www-lgge.ujf-grenoble. fr/ServiceObs.
A permanent Automatic Weather Station (AWS) set on the moraine (2700 m a.s.l., see location in Fig. 2) has recorded data since 2005. Measurements of air temperature, pressure, and relative humidity are used to compute the solar irradiance as described in Sect. 4.

Imaging radiometer system
Two digital cameras were set up next to a hut, above the glacier tongue, nearby the permanent AWS (note the location of the hut in Fig. 2). They imaged automatically the whole glacier several times a day over the period from June to Table 1. Wavelengths of maximum spectral sensivity λ c of the two cameras. The measurements of camera spectral sensitivity (Demircan et al., 2000) have been conducted using a monochromatic source and a spectro-gonioradiometer (Brissaud et al., 2004). In this study, the spectral sensitivity of each channel is approximated by a Dirac delta distribution δ λ c .

Camera
Channel 1  September. The first camera is a standard off-the-shelf camera based on a body Canon EOS 400D (22.2×14.8 CMOS sensor with 10.1 million effective pixels). The second (modified EOS 400D) is also based on a EOS 400D body, but whose sensor was modified to capture near-infrared radiation (NIR). Table 1 provides the wavelengths of maximum spectral sensitivity (Demircan et al., 2000) for each channel (RVB) of both cameras. In this study, the spectral sensitivity of each channel is approximated by a Dirac delta distribution centered in λ c . Both cameras were equipped with a Canon EF 18 mm lens to provide a field of view encompassing the whole glacier. A radiometric calibration of both cameras was conducted to determine the response function of the sensors. The latter relates the digital number of a pixel in each channel to the spectral radiance received by the camera sensor, I photo (Grossberg and Nayar, 2004). It was measured under precisely controlled illumination using Lambertian surfaces of known reflectance (Labsphere ® samples). More details about this calibration are given in Dumont et al. (2009) The cameras take photographs from the hut. Etendard peak is also indicated.

Measurements of reference albedo
In order to assess the accuracy of the maps of albedo derived from the photographs, measurements of the albedo were collected. Additional AWS were installed on the glacier during the ablation season (i.e., summer). From June to September 2008, two stations were set in the ablation and accumulation areas, respectively. During summer 2009, a single AWS was set in the ablation zone. The locations of the temporary stations are indicated in Fig. 2. These temporary stations were equipped with a net radiometer Kipp and Zonen CNR1. The two pyranometers facing upward and downward, respectively, were used to measure the incident and reflected shortwave radiation (305-2800 nm) with a 15 min time interval. The ratio of both quantities allowed the albedo of the surface to be estimated. Considering a possible tilt of the instrument with respect to the surface of the glacier and the intrinsic accuracy of the sensor in both directions (±3 % in terms of Root Mean Square Error (RMSE), Six et al., 2009), the expected accuracy of the measured albedo value is ±10 % (RMSE) (Kipp and Zonen, 2009). The albedo reference data selected for estimation of the method accuracy have an illumination zenith angle varying from 25 to 65 • .

Method
This section describes the method implemented to retrieve the albedo from two digital photographs (i.e., visible and NIR) of the glacier. In this study only clear-sky photographs were considered so as to enable a simple model to be used towards the computation of the solar irradiance at the ground level and to avoid cloud shadows on the surface of the glacier.

Preliminary processing of the digital photographs
The digital photographs of the glacier were saved in raw format to avoid unwanted distortions of the DN values measured by the camera sensor. Using the IRIS image processing software (http://astrosurf.com/buil/), the 3 channels from each photograph (i.e., visible and NIR) were isolated. Thus, measurements of spectral radiance for each central wavelength λ c (as indicated in Table 1) could be obtained.
Both cameras are placed in a box. When replacing the memory card and thus removing the cameras from their boxes, a slight tilt and/or translation may occur. Therefore, each photograph was rectified using the method described in Corripio (2004) and a Digital Elevation Model (DEM) of the glacier. To provide perfect co-registration of the visible and near IR photographs, a co-georeferencing of both photographs was conducted using manually identified ground control points and then accounting from the tilt and/or translation of the cameras.
The DEM of the glacier used in this study was generated from airborne altimetric laser measurements performed on 6 August 2008. The horizontal and vertical accuracy of the DEM is better than 50 cm. The original 1 m spatial sampling was resampled to create a DEM with 10 m pixel size using a bilinear interpolation.

Calculation of the incident radiation
The photographic measurements are related to reflected radiation only. Since the albedo is the ratio between the reflected and the incident radiant flux, it is necessary to estimate the spectral irradiance received by each pixel of the glacier. The latter was estimated using the simple parametric transmittance model described in Sirguey et al. (2009). The ground spectral incident irradiance E g (λ) is the sum of three contributions: where E s (λ) is the direct solar irradiance, E * d (λ) is the total diffuse irradiance on slope, E t (λ) is the irradiance reflected by the adjacent slopes, andθ s is the illumination angle defined as the relative angle between the normal to the surface and the solar direction (see Fig. 1b and Appendix A in this paper and Eq. (6) in Sirguey et al., 2009).

Spectral solar irradiance
is computed as the extra-terrestrial solar irradiance weighted by the sun-to-ground atmospheric transmittance T s (λ) for each pixel of the DEM. Both are calculated using a modified version of the SPCTRAL2 radiative model (Bird and Riordan, 1986) described in Sirguey et al. (2009). The total water vapour column is derived from Prata (1996) using air temperature, pressure, and relative humidity measured near the terminus of the glacier by the permanent AWS (see Sect. 3). The aerosol optical depth has been measured several times on the glacier using a sun-photometer and is assumed to be constant in the study since the measured optical depth does not vary significantly between the different days (measurements were performed at the same hour). The ozone total columns are obtained from the Ozone and Water Vapour Group (OZWV, http://www.esrl.noaa.gov/gmd/ ozwv) and this punctual measurement is considered to be valid for the whole extent of the study area.
In terms of accuracy, Sirguey et al. (2009) demonstrated that this modified version of SPCTRAL2 was in very good agreement with the more rigorous radiative model 6S (Vermote et al., 1997). Since spectral measurement of solar irradiance were not available at the study site, the solar irradiance estimated from SPCTRAL2 was integrated over the spectral domain (305-2800 nm) and compared to the measurements from the CNR1 sensor installed on the moraine. Individual measurements obtained over 10 days at a 60 min sampling interval demonstrated that the discrepancies between the two estimates was less than the inherent accuracy of the CNR1 sensor for incident irradiance (i.e., ±3 % RMSE, Six et al., 2009).

Accounting for topographic effects
As described in Eq. (5), both the diffuse environmental irradiance and the irradiance reflected by the adjacent slopes contribute to the total irradiance received by a pixel of the glacier. Both contributions involve multiple reflections between the surface and the sky, and between adjacent slopes. These effects can be addressed by means of an iterative method to estimate the incident ground irradiance . At each iteration, the ground reflectance α(λ), and ground irradiance are calculated for each pixel of the photo (see also Sect. 4.3). They are used to update the incident irradiance. The convergence of the method is discussed in Sirguey et al. (2009). Considering Eqs. (11), (23)-(25) in Sirguey et al. (2009), the ground incident radiation at each iteration (i) can be written as where the wavelength λ has been omitted for clarity. V d and C t are the sky-view and terrain configuration factors, respectively. Both were obtained from an analysis of the DEM as described in Sirguey et al. (2009). θ s is the solar zenith angle. E r and E a are the downward Rayleigh and aerosols scattering irradiance estimated using the modified SPCTRAL2 model. α e (λ) is the environmental reflectance averaged on the whole glacier.ᾱ t (λ) is the mean terrain reflectance computed in the neighborhood of each pixel by means of a 1 km moving window average. The effects of shadow were ignored and only parts of the glacier in direct sunlight were taken into account.

From measured radiance to spectral albedo
We are now able to compute the irradiance received by each pixel of the glacier. To retrieve the albedo, we must also estimate the total reflected irradiance by each pixel from the measurement made by the cameras.

Radiance measured by the camera sensor
To retrieve the albedo, it is necessary to estimate the total hemispherical radiant exitance reflected by the considered pixel from the measurements obtained from the photographs. First, the camera instantaneous field of view (IFOV) for a pixel is not infinitesimal so that the quantity measured by the camera is not truly a directional quantity. Nevertheless, the IFOV is small (less than 10 −2 sr). The effect of this solid angle was accounted for while measuring the response function (Dumont et al., 2009). Therefore, the measurement can be assumed directional. Second, the reflected radiance is attenuated by atmospheric content (mainly water vapour and aerosols) before it reaches the camera. The pixel-to-sensor transmittance was accounted for following the method described in Corripio (2004). The horizontal transmittance was www.the-cryosphere.net/5/759/2011/ The Cryosphere, 5, 759-771, 2011 computed using MODTRAN (Berk et al., 1989) with rural aerosol profile as a function of the distance between the camera and the considered pixel.
In the channel λ c , the radiance originating from a pixel and measured by the camera depends on their relative orientation. This radiance L photo (θ v ,φ v ,λ c ) is the sum of the reflections of incident radiance L g (θ i ,φ i ,λ c ) originating from all the directions of the overlying hemisphere: As explained in Sect. 4.2, the ground incident irradiance accounts for a direct and a diffuse contribution. E diff (λ) can be inferred from Eq. (5) is formed by the diffuse atmospheric radiation on slope E * d (λ) , and the terrain reflected irradiance E t (λ).
The diffuse atmospheric contribution can be considered to be isotropic (Schaepman-Strub et al., 2006) and, in a first approximation, E t (λ) is assumed to be isotropic as well. Consequently, the incident radiance originating from a direction (θ i ,φ i ) can be written as where δ 2 Consequently,

Snow and ice anisotropy
The first step to perform towards the estimation of the albedo is to retrieve the spectral albedo as defined in Eq. (2).
R(θ v ,θ s ,φ v ,λ c ) was estimated using measurements of the BRDF obtained over snow and ice with a spectrogonioradiometer as explained in Dumont et al. (2010).

From spectral reflectance to albedo
The previous section demonstrated how the spectral albedo was estimated for each pixel of the glacier in the six wavelengths corresponding to the three channels of each camera. The final objective remains to estimate the albedo (Eq. 4) for each pixel. Towards this goal, a narrowband to broadband conversion must be considered whereby the discrete values of spectral albedo are used to estimate the albedo with respect to the whole spectrum.

Database of spectral albedo for snow and ice
To perform this conversion, a database of spectral albedo for snow and ice was created using the radiative transfer code DISORT (Stamnes et al., 1988). This database was created with the Surface Specific Area (Legagneux et al., 2002) ranging from 30 to 150 m 2 kg −1 for snow while it was set to 1 m 2 kg −1 for ice. Two kinds of impurities were taken into account: mineral dust and soot. The imaginary part of the spectral refraction indices of impurities were obtained from the Gestion et Etudes des Informations Spectroscopiques Atmospheriques (GEISA) database (http://ether.ipsl.jussieu.fr/ etherTypo/?id=950&L=0). Impurity content ranging from 0 to 200 ppm was used for mineral dust and from 0 to 100 ppm for soot. The database totalized the simulated spectral albedo of 754 snow targets and 19 ice targets for 16 incident zenith angles.

Optimisation
For each pixel, the algorithm searches the database to find the spectral albedo which minimizes the distance to the six spectral albedos points estimated from the images. For most types of snow and ice, we assume that α(θ,λ)≈h(λ)g(θ,λ) where g(θ,λ) is independent of the grain size (hypothesis 1). Figure 3 shows that hypothesis (1) is acceptable up to λ≈800 nm and introduces a maximum uncertainty in terms of spectral albedo of ±5 %. Since λ max c =830 nm, this hypothesis is applicable here.
Let λ M(λ)dλ be the numerator of Eq. (4). Using Eq. (8), the radiant exitance M(λ) can be expressed as Using hypothesis (1), it comes that was calculated for all the targets in the database. For each pixel, h photo (λ i ) was computed based on the estimation of the spectral albedo obtained from the photo (Eq. 11) and the hypothesis (1). Then, the database was searched to determine the target that minimized the spectral dis- The function h DISORT of this target was used to estimate M(λ) with Eq. (13). Finally, the albedo A could be obtained using Eq. (4) for each pixel of the glacier.

Results
The retrieval algorithm was tested on the Saint Sorlin Glacier (Sect. 3) during the summers 2008 and 2009. The weather conditions enabled 12 pairs of photographs to be captured during the summer 2008 (the acquisition time was between 12:00 and 13:00 LT), while 62 pairs were acquired during the summer 2009. During summer 2009, acquisition times were set at 09:00, 12:00 and 15:00 LT in order to capture more comparison points, as well as giving an opportunity to investigate the diurnal variation and the impact of the variations of the solar incident zenith angle on the retrieval method. Reference values of albedo were measured simultaneously with the AWS described in Sect. 3. Figure 4 shows two photographs of the surface of the glacier captured by the unmodified camera on 12 June and 4 August 2009 along with the corresponding map of albedo estimated using the method described above. On 12 June, the whole surface was covered with snow. The albedo ranged from 0.5 (dirty old snow) near the terminus to 0.8 (fresh snow) at the top of the glacier. On 4 August, most parts of the surface exhibited dirty ice surfaces with albedo values as low as 0.3. Only a fraction of the surface at high elevation in the accu-mulation zone was still covered by snow whose albedo was estimated to reach 0.7.

Comparison with measurements
The scatter plot in Fig. 5 shows the values of albedo retrieved from the terrestrial photography at the locations corresponding to the AWS accumulation,2008 and AWS ablation,2008 (see Fig. 2). Figure 5 shows the albedo retrieved based on two hypothesis: (1) snow and ice were assumed to be Lambertian surfaces; (2) the BRDF effects were not ignored and the anisotropy factor described in Sect. 4.3.2 was considered. Figure 6 is similar to Fig. 5 and illustrates the comparison between the retrieved and measured albedo values during the summer 2009. During summer 2009, only one AWS was set on the glacier. Therefore, the time of the photography (i.e., 09:00, 12:00 and 15:00 LT) is illustrated in Fig. 6 with different symbols.
The summary statistics between the albedo retrieved from terrestrial photography and that measured with the AWS are shown in Table 2. The distribution of the retrieved albedo exhibits a non significant bias when the anisotropy factor is considered. The bias is larger when the surfaces are assumed to be lambertian. Similarly, the root mean square deviation (RMSE) values are generally reduced when using the correction for anisotropic effects. These statistics suggest a better accuracy of the retrieved albedo when the anisotropy factor is introduced compared to the lambertian assumption. The 99 % confidence intervals (CI) for the RMSE were computed as described in Sirguey (2009, Appendix B). When considering all sampled data (i.e., the total columns in Table 2), the CI for the RMSE values computed with and without anisotropic factor do not overlap. This demonstrates that the implementation of the anisotropy factor allowed the RMSE to be significantly reduced (at the 99 % confidence level).
In addition, Fig. 6 suggests that low albedo (ice) are over- www.the-cryosphere.net/5/759/2011/ Table 2. Mean (m) and root mean square deviation (q) of the difference between the albedo retrieved from terrestrial photography and the albedo measured at the AWS. The subscript iso corresponds to the results when ice and snow were assumed to be Lambertian surfaces. The subscript ani corresponds to the results obtained while considering the anisotropy factor. The confidence intervals on q(CI) were computed at the 99 % confidence level as explained in Sirguey (2009, Appendix B at 15:00 LT. Figure 6 also reveals five points (20 July and 20 June 09:00 and 12:00 LT, and 5 July 09:00) with very low retrieved albedo values. Finally, on 3 July 09:00 and 2 July 09:00 and 15:00 LT, the values of albedo are overestimated. Possible explanations of these features are given in Sect. 6.

Spectral distance
To evaluate the spectral consistency of the retrieval algorithm (i.e., to assess whether the retrieved spectral distribution of α photo (λ i ) is close to the simulated distribution α DISORT (λ i )), a distance d between the spectral albedo retrieved from photography and that of the corresponding surface obtained from the DISORT simulation is defined. It is formulated as the root mean square spectral difference between both quantities as follows: d ranges from 0 to 1 and is homogeneous to a spectral albedo, thus facilitating its interpretation. d=0 indicates a perfect match between the simulated and estimated spectral albedos, while d=1 indicates the largest inconsistency. For summer 2009, the average distance over all 63 measurements was d=4.90×10 −2 (d=4.65×10 −2 for summer 2008 on 12 photographs).

Discussion
This section provides explanations of the results presented in Sect. 5 and discussions on the accuracy and sensitivity of the retrieval method.

Results interpretation
The results presented in Table 2, Figs. 5 and 6 illustrate the performance of the retrieval method and demonstrate that the implementation of an anisotropic model of reflectance improves the accuracy of the estimated value of albedo. In Fig. 6, three data points collected on 3 July 09:00 LT and 2 July 09:00 and 15:00 LT exhibit an overestimated albedo. It is believed that an error while georeferencing the photographs could have caused this difference. The glacier was visited on 3 July and the AWS was located close to the transient snow line (Fig.7). The latter line is the boundary between remaining snow and ice and is close to the equilibrium line at the end of the ablation season (line of zero mass balance). Thus, the difference in albedo above and below the AWS was not negligible. Consequently, a slight error while georeferencing the photographs could have resulted in the retrieval of the albedo of old snow (instead of dirty ice) at the pixel corresponding to the location of the AWS. This was supported by the strong albedo gradient found at this location on the map of albedo. In addition, the operator noticed the presence of a large amount of liquid water at this location which might also have affected the accuracy of the algorithm.
The presence of high clouds (cirrus) could explain the underestimation of the albedo on other dates (e.g., 20 July 09:00 and 12:00 LT, 20 June 09:00 and 12:00 LT and 5 July 09:00 LT, see Fig. 6). High clouds reduce the incident irradiance at the ground level although they are ignored by the atmospheric transmittance model. Thus, E g could have been overestimated and α underestimated. To test this hypothesis, the level of cirrus reflectance given by MODIS level 2 -Cloud Products (1 km spatial resolution) was investigated. From all the images that matched the photographs (day and time), only one day was found to exhibit a high value of cirrus reflectance (20 June 2009). On 5 July 09:00 LT, clouds were visible in the NIR photography. On 20 July no high cloud was visible on satellite data. Nevertheless, snow fell on 18 July. On 19 July the algorithm provided suitable estimation since the glacier was covered with fresh snow. However, on 20 July patches of melting snow and ice were present on the lower part of the glacier.  of albedo revealed a high interpixel variability in the direct neighborhood of the AWS (α is ranging from 0.4 to 0.6 from pixel to pixel). At 12:00 LT, the interpixel variability reduced (α ranging from 0.4 to 0.5). However, it is believed that the CNR1 (10 m-diameter circular detection at the surface) was still facing melting snow (the measured value of albedo was 0.572) although the pixel was mostly formed by bare ice. This could explain the discrepancy between both values of the albedo.
6.2 Sources of error

Georefererencing
As noticed above, the rectification of the photographs can introduce errors in terms of pixel location. The maximum georeferecing error was found to be less than ±50 m. In addition, the georeferencing error is sensitive to the accuracy of the Digital Elevation Model. Elevation and macroscale roughness of the surface of the glacier (often created by channels of   Table 3. For both scenarios, the difference in terms of retrieved albedo is smaller than the uncertainties introduced by the method. The differences are larger in June, when the albedo is larger (snow is present instead of ice), than in August. It appears that the retrieved albedo value is not very sensitive to reasonable errors in the DEM (i.e., errors typically less than 4 m in elevation).

Atmospheric modelling
The modelling of the atmospheric irradiance in the modified SPCTRAL2 model is another source of uncertainty. Under clear-sky, this error was evaluated by Sirguey et al. (2009) to be less than ±10 % of the total irradiance. Obviously, the error can increase in the presence of high clouds as noted above.
The other underlying hypothesis in the method is the fact that diffuse irradiance, atmospheric as well as terrain reflected, is considered isotropic. However these two diffuse radiances originate from exclusive regions of the overlying hemisphere: the region where the sky is visible for the atmospheric diffuse contribution and the region where the terrain is visible for the terrain reflected contribution. Consequently, a more rigorous approach would be to keep both components separated and to modulate them with the sky view factor. Nevertheless, we believe that the isotropic assumption is justifiable with respect to the limited amount of radiation that is involved compared to the direct radiance provided that the ilumination angle is not too high (typically less than 60 • ), which is the case on most of the photographs.

Anisotropy
Uncertainties can also originate from the correction of the anisotropy. It is true that the correction applied by means of a simple anisoptropic factor is relatively rough. Full BRDF models (Tedesco and Kokhanovsky, 2007;Lyapustin et al., 2009) could have also been used but they require more information about the targets. Nevertheless, as the surface roughness and elevation of the glacier is changing over the summer, an exact knowledge of the viewing and illumination angles is not achievable. Therefore, the implementation of a more accurate BRDF correction is believed unnecessary as the roughness of the surface and the variations of elevation are not accurately known.

Spectral content
Finally, it should be stressed that no spectral information was measured for wavelengths greater than 830 nm (Table 1). Additional uncertainties may result from the fact that the effect of grain size and shape on spectral albedo become significant at wavelengths greater than 1000 nm (Warren, 1982). To evaluate this uncertainty, two spectral albedos were selected in the DISORT database of snow (see Sect. 4.4.1) with same content of impurity but different values of Specific Surface Area. The values of albedo were calculated for maximum and minimum SSA using typical ground irradiance on Saint Sorlin Glacier during summer. Maximum and minimum SSA were set so that |α(SSA max )−α(SSA min )|≈0.05 albedo unit at 800 nm (i.e., d=0.05, see Sect. 5.3). The maximum relative difference in terms of albedo was found to be less than ±2 %. Thus, it is argued that the fact that no spectral information is available for λ>830 nm does not imply an error larger than ±2 % on the final retrieved albedo value.

Conclusions
This study presents a new method to retrieve surface glacier albedo from terrestrial photographs while taking into account the anisotropic reflectance of snow and ice. Narrowto-broadband conversion is performed using the DISORT radiative model (Stamnes et al., 1988). The accuracy of the method is estimated to be better than ±10 % under clearsky conditions. However, the presence of high clouds and high intra/inter-pixel variability negatively affect the accuracy. The method described in this study allows retrieval of broadband albedo under clear sky conditions. Nevertheless, albedo varies significantly with cloud cover (Dozier et al., 2009). Thus a future refinement of the method will be to implement a new solar irradiance model which takes into account the effect of the cloud cover and allow retrieval of the albedo under cloudy conditions. Finally, the georeferencing process is one of the main sources of uncertainty.
The range of albedo reference measurements is quite large (25 to 65 • ). Thus, the method can also be easily applied to www.the-cryosphere.net/5/759/2011/ The Cryosphere, 5, 759-771, 2011 polar glaciers. Probably the anisotropy correction at larger zenith angle and the surface roughness will be crucial. For illumination zenith angles larger than 65 • , the anisotropy factor used in this study might need further improvements. This method allows the spatial and temporal variations of surface albedo to be monitored closely. This has the desirable potential to contribute towards better modeling of glacier mass balance. Indeed the albedo maps retrieved from the photographs can be assimilated into an energy balance model to accurately simulate the spatial distribution of the glacier mass balance during the ablation season .
The method is also applicable to satellite data such as MODIS. Indeed, currently available MODIS products have often large spatial and temporal resolution and are validated on rather flat areas (Dozier et al., 2009). The method presented in this study applied to MODIS radiances, thanks to accurate topography and anisotropy corrections, will probably allow retrieval of broadband albedo on highly rugged topography at 250 m resolution. Although satellite data have often poorer spatial resolution than terrestrial photography, their synoptic view provides a larger spatial coverage. The application of the method to satellite data will allow monitoring of several glaciers, while avoiding the setup and maintenance of specialized digital cameras.

Angles definitions
Let (S) be the frame of reference attached to the glacier where the z-axis is vertical and x-axis is northward. Let (S) be the frame of reference attached to a pixel of the glacier where the z-axis is the normal to the pixel, n, and the x-axis is aligned with the azimuth of the sun.
Let s be the vector defining the direction of the sun, v be the vector defining the direction of the observation. Tilded angles are defined in the reference frame (S).
The illumination angle,θ s , is the relative angle between s and n so that cosθ s = s × n = cosθ s cosθ n + sinθ s sinθ n cos(φ s − φ n ). (A1) Similarly, the relative zenith observation angle,θ v , is defined by cosθ v = v × n = cosθ v cosθ n + sinθ v sinθ n cos(φ v − φ n ). (A2) Finally the relative azimuth,φ v , between the sun and the observation can be computed using projections of s and v (s andṽ) on the surface tangent to the pixel given bỹ v = v − (v × n)n ands = s − (s × n)n. Consequently, cosφ v =ṽ ×s ṽ s (A4) which corresponds to cosφ v (A5) = cosθ v cosθ s +sinθ v sinθ s cos(φ v − φ s )−cosθ v cosθ s sinθ s sinθ v .