High-accuracy measurements of snow Bidirectional Reflectance Distribution Function at visible and NIR wavelengths - comparison with modelling results

Abstract. High-accuracy measurements of snow Bidirectional Reflectance Distribution Function (BRDF) were performed for four natural snow samples with a spectrogonio-radiometer in the 500–2600 nm wavelength range. These measurements are one of the first sets of direct snow BRDF values over a wide range of lighting and viewing geometry. They were compared to BRDF calculated with two optical models. Variations of the snow anisotropy factor with lighting geometry, wavelength and snow physical properties were investigated. Results show that at wavelengths with small penetration depth, scattering mainly occurs in the very top layers and the anisotropy factor is controlled by the phase function. In this condition, forward scattering peak or double scattering peak is observed. In contrast at shorter wavelengths, the penetration of the radiation is much deeper and the number of scattering events increases. The anisotropy factor is thus nearly constant and decreases at grazing observation angles. The whole dataset is available on demand from the corresponding author.


Introduction
Snow covered areas on earth reflect more solar radiation than any other surfaces.Snow albedo is an important parameter to accurately compute the radiation budget of regions covered by seasonal or permanent snow and thus has a significant influence on the earth radiation budget (Jin et al., 2008).
An accurate retrieval of this variable and a better understanding of its variations with zenith angle, wavelength and snow physical parameters is required to study the influence of snow cover changes on climate (Xie et al., 2006).
Considering the high spatial and temporal variability of albedo and the fact that most snow covered areas are difficult places to reach to perform field measurements, remote sensing is the most suitable tool to determine spatial and temporal variability of snow albedo.Nevertheless, most of remote sensing sensors measure the reflected radiation into a few wavelength bands and a particular angle (bi-conical reflectance, Nicodemus et al., 1977) instead of the albedo.Since snow is not a Lambertian reflector (Lyapustin et al., 2009;Li et al., 2007;Warren, 1982), the conversion from bi-conical reflectance, so measured by satellites, to spectral albedo is not straigthforward.To convert the bi-conical reflectance into a spectral albedo useful for radiation budget calculation, the angular distribution of the reflected radiation on snow has to be known.In lighting and viewing solid angles are infinitesimal, the biconical reflectance is named bidirectional reflectance and its angular distribution over the upper hemisphere the Bidirectional Reflectance Distribution Function (BRDF) (Nicodemus et al., 1977).This angular distribution is described by the anisotropy factor, the BRDF normalized by the spectral albedo (Warren et al., 1998;Li et al., 2007).
The intent of this study is to present highly accurate measurements of snow BRDF performed in a cold room using a spectrogonio-radiometer and to compare them with previous field measurements.Definitions related to albedo and reflectance are given in Sect. 2. Previous field measurements and modelling studies are quoted in Sect.3. Four different types of snow are sampled in order to study the sensivity of BRDF to snow grain size, shape and impurity content.A wide range of lighting and viewing geometries are explored.These measurements are one of the first sets of direct, i.e. not under natural illumination, snow BRDF values over a large spectral range.They are also one of the first investigations over a large range of incident lighting configurations with incident zenith angle varying from nadir to 60 • (Sects.4, 5 and 6).We compare these measurements with modelled BRDF in order to investigate the accuracy of radiative transfer models and to understand the scattering phenomena observed in the measurements.In this study two different radiative transfer models are applied: the SnowRAT (photon tracing) model (Picard et al., 2008) and Mishchenko model (Mishchenko et al., 1999) (Sect.7).Finally, Sect.8 provides a discussion and Sect.9 conclusions.
This study considers flat surfaces and snow grains are assumed to be randomly oriented.Hence the BDRF is symmetric with respect to the principal plane which contains the incident beam and the normal to the surface, BRDF only depends on the relative azimuth, φ=|φ i −φ v |.Consequently in the following the incident beam azimuth angle is 0 and φ∈[90 • ; 270 • ] defines the forward part of the hemisphere.
The spectral albedo α(θ i ,λ), or directional-hemispherical reflectance (Nicodemus et al., 1977), is the ratio of the reflected radiation in the whole upper hemisphere over the incident collimated irradiance at a given wavelength.
In order to isolate the angular dependence of snow surface reflection, the anisotropy factor (Nicodemus et al., 1977), R(θ i ,φ,θ v ,λ) is constant and equal to unity in the case of a perfectly Lambertian surface.
In most field studies, BRDF values are not accessible by measurement since natural light is not collimated but includes both direct solar and diffuse radiations.The Hemispherical Directional Reflectance Function (HDRF) h, is the accessible parameter, which means that the incident beam is integrated over the whole incident hemisphere and includes both direct and diffuse irradiances.Nicodemus et al. (1977).
3 State of the art

Field studies
Many field studies of snow reflectance properties have been performed.To our knowledge, only a limited number of direct snow BRDF measurements are available since natural Atmos.Chem.Phys., 10, 2507-2520, 2010 www.atmos-chem-phys.net/10/2507/2010/light is not collimated.Indeed most of the studies referenced below give access to HDRF and not BRDF.The latter allows to calculate reflected radiance for any given incident sky radiance distribution while HDRF is applicable to the specific conditions of illumination during the observation.Nevertheless HDRF is close to BRDF in the infrared part of the spectrum under natural conditions (Li and Zhou, 2004).Warren (1982) gave an overview of previous measurements of snow reflectance.Leroux et al. (1998) and Sergent et al. (1998) investigated the influence of snow grain shape and size on HDRF in both cold laboratory and field.Warren et al. (1998) and Grenfell et al. (1994) studied the effect of macroscale roughness on snow HDRF using measurements at South Pole Station at 3 different wavelengths.They concluded that macroscale roughness significantly influences HDRF patterns.Aoki et al. (2000) measured in the field HDRF from 350 to 2500 nm and analysed the effects of grain size and impurities on it.Several studies presented different field measurements of HDRF over several parts of the solar spectrum: [350-2500] nm (Painter andDozier, 2004), [350-1050] nm (Li and Zhou, 2004;Bourgeois et al., 2006) and[390-1070] nm (Peltoniemi et al., 2005) and studied the influence of solar zenith angle and various snow properties (grain size and shape, wetness, impurity, depth, density).Hudson et al. (2006) described measurements of HDRF at Dome C, Antarctica at solar zenith angles from 51 • to 87 • .Kaasalainen et al. (2006) performed accurate measurements of the backscattering peak on numerous snow samples.These field studies illustrate the main patterns of the snow anisotropy factor and of its variations with the physical properties of snow and with viewing and lighting angles, but remain limited by the accessible zenith angles and uncollimated incident radiation under natural conditions.Wiscombe and Warren (1980) first introduced a model for the computation of spectral albedo based on Mie theory and the δ-eddington method.Leroux et al. (1998Leroux et al. ( , 1999) ) used adding-doubling method to calculate snow HDRF and spectral albedo.Warren et al. (1998) underlined the importance of an accurately modelled phase function for the computation of snow BRDF.Painter and Dozier (2004) and Li and Zhou (2004) used respectively DIScrete Ordinates Radiative Transfer (DISORT) (Stamnes et al., 1988) and addingdoubling method with equivalent spheres of equal volumeto-surface-area ratio as snow grains.They noticed that a precise computation of single scattering parameters (i.e.single scattering albedo and the phase function) is essential to simulate accurate BRDF.Aoki et al. (2000) compared Mie theory and Henyey-Greenstein semi-empirical phase function to model single-scattering parameters and concluded that spectral albedo can be accurately simulated using equivalent spheres whereas BRDF cannot.Grenfell and Warren (1999) showed that simulating one non-spherical particle by a collection of independent spheres of same total volume-tosurface-area ratio leads to accurate retrieval of single scattering parameters.Mishchenko et al. (1999) presented a model for the computation of snow BRDF based on an analytic solution of the radiative transfer equation and an approximation of the phase function.This model is applicable for any shapes of particles.Kokhanovsky and Zege (2004) presented an asymptotic solution of radiative transfer theory adapted to snow and able to deal with fractal and spherical particles.They concluded that fractal particles are more appropriated for simulating snow BRDF than spheres.Kokhanovsky et al. (2005) compared the results of their asymptotic model with in situ measurements and concluded that the accuracy is reduced in the principal plane and at high observation angles.Xie et al. (2006) compared three different radiative transfer models (DISORT, adding-doubling and Mishchenko model) and two truncation methods of the forward peak (δ-eddington and δ-fit).They concluded that only an accurate computation of single scattering albedo, ratio of scattering efficiency to total light extinction i.e. scattering and absorption, is essential to account for the influence of grain size on BRDF.Besides, in order to account for the influence of grain shape, both single scattering albedo and phase function should be accurately simulated.Picard et al. (2008) used a photon tracing model to compute snow albedo for several grain shapes in the near IR.Additionally Jin et al. (2008) used a coupled snow-atmosphere model to generate anisotropy factor and spectral albedo for layered snowpack and validated their approach with measurements of Hudson et al. (2006).All these studies lead to the conclusion that spectral albedo can be accurately modelled but that theoreticals difficulties linked with the non sphericity of snow grains still remain and limit accurate modelling of snow BRDF.

Spectrogonio-radiometer
The BRDF has been measured using the spectrogonioradiometer developed at the Laboratoire de Planétologie de Grenoble, France.A comprehensive description of the device is given by Brissaud et al. (2004) and Bonnefoy (2001).The sample is illuminated by a monochromatic light with a spectral width from 0.2 nm to 0.6 nm.The incident zenith angle varies from 0 to 80 • with a beam resolution of ±0.1 • .The viewing zenith angle varies in the same range and the azimuth angle takes any value from 0 to 180 • .
At nadir incidence, the illumination pattern at the sample surface is circular, with a 200 mm diameter.The spatial variations of the light intensity inside the area viewed by the detector are typically less than 1% at nadir.The detectors have an half angle field-of-view of 2.05 • and a circular observation pattern of 20 mm diameter at nadir, but a larger elliptic pattern at other incident angles.
As noticed in Schaepman-Strub et al. (2006), BRDF as defined in Eq. 1 cannot be directly measured since it requires an infinitesimal solid angle of observation.Thus the quantity measured by the spectrogonio-radiometer is the directionalconical reflectance (Schaepman-Strub et al., 2006).Nevertheless, since the detectors fields of view are small we consider in the following that our measurement are very close to BRDF.
The device is located in a cold room at −10 • C to allow measurements on snow.

Snow samples
Four samples of various types of snow (S1, S2, S3, S4) have been collected at various locations in the French Alps in January 2008.The samples are cylindrical (30 cm diameter, 12 cm deep) and large enough to minimize edge effects even at visible wavelengths.Test experiments have been conducted at 630 nm on a transparent cubic sample holder (29.5×29.5×16.5 cm 3 ) filled with artificial snow.The results show that side losses are less than 0.5% and bottom losses less than 1% at 0 • incident zenith angle.Consequently photon losses in the 2 cm observation pattern are estimated to be less than 0.1%.The reflecting sides of the sample holder further decrease these losses.
The samples, except S4 dedicated to test the temporal evolution of snow during measurements, were stored at −10 • C during at least one week before being measured in order to allow thermal stabilization and to avoid metamorphism during the measurements.During the first hours in the cold room, the wet snow samples (S2, S3 and S4) refroze.
Digital pictures have been taken in order to characterize the grain shape and size for each sample.No impurities content measurements were performed and only S2 contains a high quantity of impurities visible by eye.Actually, S2 was taken near Argentière village (Mont-Blanc valley).Table 1 shows the sample characteristics.Grains digital photographies show that r(S1)≈r( S3)<r(S4)<r(S2) where r represents the effective radius of the sample.

BRDF measurements
For each sample (except S4), a complete set of radiance measurements has been performed as follow.The spectral range covered is 500 to 2600 nm with a 20 nm step.
Three different incident angles, θ i , have been chosen (0 • , 30 • and 60 • ) in order to study the effect of lighting zenith angle on BRDF.For each incident angle, the viewing zenith angle, θ v , takes different values: 0 • , 30 • , 60 • and 70 • and the relative azimuth, φ, is 0 • , 45  10, 2507-2520, 2010 www.atmos-chem-phys.net/10/2507/2010/Measurements at larger incident or observation zenith angles were not performed since edge effects due to the size of the sample holder are too significant for these configurations.
Radiance measurements on S4 were limited to a single geometry but repeated for 24 h to investigate the stability of the measurements.

Reference measurements
To convert spectrogonio-radiometer measurements of the reflected flux into reflectance values, we divide the snow measurements by the flux reflected by reference surfaces for which spectral albedo and BRDF are known.For visible and near-IR wavelengths, the reference surface is a Spectralon ® , a nearly perfect Lambertian reflector.For IR wavelengths, longer than 2440 nm, an infragold ® sample is taken as reference since Spectralon ® is unsuitable at these wavelengths.The relative accuracy of the reflectance measurements is better than 1% using fully calibrated references (Bonnefoy, 2001).

Shadow and geometric limitations
Due to the size of the detectors, no measurement can be performed in the backscattering direction (i.e. in the nearly same direction as the incident beam, (θ i ≈θ v , φ≈0)).
In order to convert the BRDF measurements into spectral albedo and anisotropy factor, we assume that the BRDF is symmetric with respect to the principal plane (Hudson et al., 2006) (for azimuths from 180 • to 360 • ) and perform a linear interpolation in cos(θ v ) and φ for our measurements over the whole observation hemisphere.Extrapolation of the measurements has also been performed firstly from 70 • to 90 • observation zenith angle for all the incident zenith angle and secondly to fill the blanks due to shadows of detectors.
Results of simulation with Mishchenko's model (Mishchenko et al., 1999) and SnowRat model (Picard et al., 2008) have shown that the uncertainties on the spectral albedo values resulting from interpolation and extrapolation is less than 2% except at very high absorption value i.e. spectral albedo smaller than 0.01.An estimation of the resulting uncertainties has been plotted in Fig. 2. Spectral albedo values used in R charts are indicated in the legend of each figure.

Temporal evolution of the sample during measurements
The measurements for one sample lasted approximately 33 h.
To check that the snow structure did not change during the acquisition we performed reflectance measurements on fresh snow immediately after being collected (S4) and we repeated the same measurements 24 h later.Comparing both Fig. 2. S3 spectral albedo α(λ,θ i ) calculated from BRDF measurements for three incident zenith angles.Vertical lines are located at spectral albedo minima.Gray areas represent an estimation of the uncertainties of albedo values due to interpolation and extrapolation over the measurements of the whole hemisphere.The uncertainties have been evaluated using SnowRAT (Picard et al., 2008) and Mishchenko's models (Mishchenko et al., 1999).
reflectances (not presented here), we can notice that the absolute difference in reflectance is smaller than 0.01 for wavelengths shorter than 1 µm.For longer wavelengths, the absolute difference is larger and reaches 0.025 at 1.4 µm.This later wavelength corresponds to an absorption secondary mimimum where the reflectance sensitivity to grain type is maximum (Wiscombe and Warren, 1980).
As a conclusion and since fresh snow is more subject to metamorphism than aged snow, metamorphism has only a moderate effect during the 33 h of our measurements.

General patterns of snow spectral albedo, BRDF and anisotropy factor
Only S3 measurements are presented in this section as they are representative for the other samples.Spectral albedo, α(λ, θ i ), for the three incident angles are plotted in Fig. 2. It takes the highest values in the visible and decreases at longer wavelengths with 4 remarkable local maxima and 4 secondary minima due to ice absorption bands (1.03, 1.26, 1.5 and 2 µm).This plot also shows that spectral albedo increases at all wavelengths with incident zenith angle.
Figure 3 shows the spectral anisotropy factor, R(λ, θ i , θ v ), as a function of wavelength, for three different observation angles, with a fixed illumination angle (θ i =30   spectral albedo.A low spectral albedo corresponds to large variations of R with observation angles.In the visible and up to 1.2 µm, R is nearly constant and close to unity.At longer wavelengths, R strongly diverges from unity. Figure 4 shows R values at two wavelengths selected for their differences in absorption.It underlines the increase of anisotropy with absorption and shows three important scattering features with respect to incident angle.(1) At nadir lighting, snow reflectance is nearly Lambertian.However, the anisotropy factor is not fully circularly (varying only with θ v ), as it should be for perfectly horizontal sample.(2) At 30 • incident angle, and 0.6 µm, R shows a forward scattering peak at (θ v =30 • , φ=180 • ).This feature is referred as darkening at grazing angles in the following since R is decreasing at limb in the forward direction.R maximum increases and shifts to larger viewing angles as wavelength increases.(3) At higher incident angle (60 • ), the forward scattering peak becomes sharper and stronger and is observable at both wavelengths.

Influence of snow physical properties on the anisotropy factor
To investigate the influence of snow physical properties (size and shape of grains and impurity content) on snow anisotropy factor, we compute the ratio of R for two different samples: R(S1) R(S3) in Fig. 5 and R(S1) R(S2) (Fig. 1 in the supplement: http://www.atmos-chem-phys.net/10/2507/2010/acp-10-2507-2010-supplement.pdf) at several wavelengths and for three incident zenith angles.
Figure 5 shows that the R ratio varies by less than 10% at wavelengths shorter than 1 µm whatever the incident angle is.However at wavelengths greater than 1 µm, both magnitude and angular patterns of R change significantly as a function of wavelength.As an example, at 1.5 µm, 30 • incident angle and (θ v =30 • , φ=180 • ), S1 anisotropy factor is 1.5 times higher than S3 anisotropy factor.At 1.5 µm, the variability of R(S1) R(S3) seems to be larger at 60 • , especially in the side and backward directions.As for the forward direction, the maximum value is larger at 30 • but the ratio migth be larger at 60 • if observations were extended at 80 • viewing angle.Moreover, concerning R(S1) R(S2) , the variability of the ratio increases with incident zenith angle and seems to indicate that for small and elongated grains, the forward scattering peak is higher in magnitude than for large rounded grains.As a conclusion, Fig. 5 shows that snow grains shape and size have little impact on the shape of the BRDF in the visible and up to 1 µm whereas for longer wavelengths the effect is much stronger.
This point is reinforced by Fig. 6 which presents R for the three snow samples at 1.5 µm (one of the absorption maxima) at 30 • incident angle.The shape and magnitude of the forward scattering peak clearly depends on grain size and shape.Indeed for S1 (broken dendritic crystals), R presents two separated maxima at (θ v =30 • , φ=180 • ) and (θ v =70 • , φ=180 • ).Only one maximum is observed for S3.For S2, we are not able to assume whether there is only one maximum or whether the angular sampling of the measurements does not allow observing two separated maxima even if the measurements at 2.5 µm in Fig. 5 in the suppplement: http://www.atmos-chem-phys.net/10/2507/2010/acp-10-2507-2010-supplement.pdfseem to validate the hypothesis of two separated maxima.However, for S3 composed of rounded grains and mixed forms, the forward scattering peak is confined at grazing observation angles.The magnitude of R maximum also varies with grain shape and size and is much greater (2.2) for S3 and S1 than for the other sample (1.8).
Figure 7 plots for each sample, as in the work of Hudson et al. (2006), the anisotropy factor in the forward scattering peak (θ v =70 • , φ=180 • ) as a function of spectral albedo at 30 • and 60 • incident angles.It shows that R and α accurately follow a power law relationship which depends on sample and illumination angle.This relationship only sligthly degrades at very low albedo values.

Comparison between measurements and modelling results
This modelling study aims at testing the consistency of the measurements and at understanding the physical processes determining the angular location and intensity variations of R maximum that appears in Figs. 4 and 6.The mean grain sizes used as inputs are set to 0.1, 0.4 and 1 mm, typical values for our samples.However we used empirical grain size distributions as our grain size measurements are too rough to be suitable as inputs of the models.
We consider two models for the calculation of snow BRDF: a model that uses an analytic solution of the radiative transfer equation (Mishchenko et al., 1999) and a photon tracing model, SnowRAT (Picard et al., 2008).SnowRAT is a discrete model whereas Mishchenko model is a continuum model.The complementarity of both models allows to understand the various phenomena observed in our measurements.

Results with Mishchenko model
Mishchenko model (Mishchenko et al., 1999) accuracy is maximum at wavelengths associated with low or intermediate values of absorption.
www.atmos-chem-phys.net/10/2507/2010/Atmos.Chem.Phys., 10, 2507-2520, 2010 To limit the computation time, we only performed simulation for spheres.Figure 8 presents the results of the simulations for a power law distribution of sphere radius with an effective radius of 100 µm (Eq.22, Mishchenko et al., 1999).General patterns of R are very similar with the measurements.Calculated R increased at longer wavelengths and reaches 2.5 at 1.5 µm in comparison with 2.2 in the measurements in Fig. 6.The darkening at grazing angles is visible at 0.6 and 1.02 µm.The effect is however stronger in the model.In Fig. 8, the decrease of R from its maximum value to (θ v =70 • , φ=180 • ) value is 0.2 at 0.6 µm and 0.1 for the same geometry at 1.02 µm.In comparison in the measurements, the corresponding decrease of R are respectively 0.10 and 0.05.
One can also notice in Fig. 8 the crescent centered around the incident direction that appears in the backward direction at all wavelengths.

Results of the photon tracing model
SnowRAT (Picard et al., 2008) ig. 7. Anisotropy factor R(θ v =70 • , φ=180 • ) versus spectral albedo (log-log plot) for the three samples and two incident angles (30 • and 60 • ).Each point on the chart belongs to a wavelength for one sample and one incident zenith angle.
variations and range of anisotropy factor, the same crescent predicted in the backward direction).The results obtained for random cylinders are presented in Fig. 9, for two radii (0.4 mm and 1 mm).Most of the main characteristics of the measurements (Fig. 4) are well reproduced i.e. the anisotropy generally increases with increasing wavelengths; darkening at grazing angles appears at 0.9, 1 and 1.3 µm; R at (θ v =70 • , φ=180 • ), at 1 µm, differs by 0.2 from its primary maximum and for the secondary maximum of absorption at 1.5 µm, two R maxima appear at (θ v =30 • , φ=180 • ) and at (θ v =70 • , φ=180 • ), as for S1 (Fig. 6a).
In summary, anisotropy factors obtained by measurements, radiative transfer and photon tracing models show strong similarities: forward scattering, darkening at grazing angles and double maxima of the anisotropy factor.However slight differences between measurements and models still exist.As an example, the photon tracing model predicts that the anisotropy for 0.4 mm radius is globally smaller than for greater radii (1 mm).The measurements (Fig. 6) present more contrasted R for small particles (S1) than for large grains (S2).

General variations of spectral albedo and anisotropy factor
Snow spectral albedo increases with incident zenith angle (Fig. 2) as explained by Warren (1982).At near nadir incident illumination, photons escape the snowpack with a lower probability than at grazing incident angle.The spectral www.atmos-chem-phys.net/10/2507/2010/Atmos.Chem.Phys., 10, 2507-2520, 2010 albedo is thus lower.This allow to distinguish two cases: in the first case (θ i ≥60 • ), photons stay near the surface and single scattering prevails with respect to multiple scattering; in the second case (θ i ≈0 • ), photons penetrate deep into the snow-pack and the number of scattering events is high before escaping or absorption.Furthermore, as noticed earlier, spectral albedo increases with incident zenith angle at all wavelengths including visible wavelengths.This quite unsual feature is most probably due to the fact that snow sample contains impurities.Thus visible wavelengths present the same pattern as more absorbing wavelentghs due to absorption caused by impurities.It largely differs from the behaviour of pure snow which is highly transparent at visible wavelengths.
Besides, in Figs. 2 and 3, snow spectral albedo and anisotropy are anti-correlated; spectral albedo is globally decreasing with wavelength and the angular distribution of R is globally more contrasted.R variations are correlated with absorption which is proportional to the imaginary part of ice refraction index presented in Warren et al. (2008).In addition anisotropy presents maxima in the absorption bands where spectral albedo presents minima.While absorption increases, the probability for a photon to be absorbed is higher.Thus, spectral albedo decreases and most photons which are reflected have only undergone a limited number of scattering events near the surface.Consequently, as absorption increases, the number of scattering events decreases.R is then mostly controlled by single scattering parameters of individual snow grains which are strong forward scatterers (Xie et al., 2006) and R increases.Furthermore, Fig. 7 corroborates the fact that R is physically related to the absorption as explained by Hudson et al. (2006).
We now propose some interpretations on the angular variations of the maximum of R as a function of wavelength and incident angle.Two phenomena are mainly observed (1) darkening at grazing angles and (2) forward scattering peak.
(1) At wavelengths shorter than 1 µm, R patterns show darkening at grazing angles in situations of near-nadir incidence (0 • , 30 • ) (Fig. 4).Darkening at grazing angles also appears on model results (Figs. 8 and 9) whatever the shape of the grains is.This effect is as well noticeable in Fig. 3a in Hudson et al. (2006) but only at non-forward directions.The absorption is small at these wavelengths and the number of scattering events that a photon undergoes before escaping or absorption is high.To understand this effect, the source function S, is a useful tool.S is the potential related to the emergent intensity along the optical path.S at any depth is the result of scattering, both upward and downward, out of an infinitesimal volume at that depth.The source function generally decreases from the surface to depth, except in the uppermost layers where it increases with depth because downward diffuse radiation approaches zero very close to the surface.Thus the source function initially increases with depth from the surface, reaches a maximum just below the surface and then steadily decreases.Consequently at grazing observation angles, the energy emerges from the region located above the maximum in the source function.Thus radiance is lower than at near vertical observation angles in which case the energy emerges from the region located deeper in the snowpack, including the maximum of the source function.This explanation is independent of the phase function which depicts the density of probability for a photon that impacts the grain to be deviated of a given scattering angle from its original direction.Measurement and model results confirm this assumption.Darkening at grazing angles is due to multiple scattering and occurs whatever are the size and shape of grains but solely for incident angles close to nadir and at wavelengths with low to moderate absorption.
Another explanation to darkening at grazing angles is proposed based on the photon concept.This effect implies that the highest probable outgoing direction for a photon that is scattered out of the snow-pack is near vertical.This is explained as follow: considering a photon, at a given depth in the snowpack, known to escape the snowpack in the future.The most probable trajectory is the one for which it will be least scattered.This trajectory is the shortest way.Consequently for a nearly isotropic photon flux inside the snowpack (multiple scattering), if it escapes at the surface, the most probable path, is a near vertical path.
(2) At wavelengths larger than 1 µm or at large incident zenith angle, R patterns show a strong forward scattering peak.This appears both in measurements (Figs. 4 and 6) and in model results (Figs. 8 and 9).At these wavelengths, the ice absorption is significant and scattering mainly occurs close to the surface.Consequently, R patterns are mostly controlled by single scattering properties of grains and especially the phase function.The phase functions for different grain shapes (spheres, random particles, columns ...) are characterized by a maximum around 0 • scattering angle (Xie et al., 2006;Warren, 1982;Mishchenko et al., 1999).This means that the most probable direction of emergence for a photon that intersects an ice particle is the straight direction.Since in case of strong absorption, a photon that escapes the snowpack has only undergone a limited number of scattering events, it will most probably escape at viewing angle close to 90 • in the forward direction when the incident angle is large enough (≥30 • ).This explains the forward scattering peak that appears at high viewing angles (≥60 • ) on R(θ v ,φ) charts.
In addition the strong correlation observed in the log-log plot between R(θ v =70 • ) and α (Fig. 7) results in an analytic relationship of the type R=Aα B , where A and B both vary with grain shape and size and incident angle.In practice, this relationship can be a useful tool to parametrize R (Hudson et al., 2006).In a more theoretical way, the meaning of the relationship should be investigated further.
As a conclusion, two cases can be distinguished: (1) low value of absorption or near nadir incident angle, photons penetrate deep into the snowpack and R is mostly controlled by multiple scattering; (2) strong absorption or/and high incident angle, photons are scattered near the surface and R is mainly determined by the particle phase function.

Effect of grain size and shape on snow anisotropy factor
Section 7.2 points that photon tracing simulations and measurements give contradictory results concerning R variations with grain size and shape.Figure 5 shows almost no variation of R between samples at wavelengths smaller than 1 µm.As explained by Painter and Dozier (2004), at these wavelengths the influence of grain size and shape is limited due to low absorption and the large number of scattering events.Consequently we believe that the slight differences observed for the ratio R(S1) R( S3) might be due to different impurity contents.At wavelengths larger than 1 µm, the ratio of measured anisotropy factors markedly differs from unity.Variations of R for S1 and S3 are higher than for S2 (the coarsest grains).The anisotropy simulated with SnowRAT (Fig. 9) is, in contrast, greater for larger grains and the forward scattering peak is stronger.The two results disagree.In addition, HDRF measurements by Painter and Dozier (2004) indicate that R is greater for large grains than for fresh snow.In Bourgeois et al. (2006) the forward scattering peak is also larger for large grains.
A careful study of the results presented by Xie et al. (2006) underlines the fact that single-scattering albedo decreases as grains size increases whatever is the shape of grains.Consequently, R becomes larger with increasing grain size as long as grain shape does not change.However, the phase functions in Fig. 5 in Xie et al. (2006) vary significantly with grain size.This change does not take the same direction depending on the shape of the grains.
At constant grain shape, R increases with grains size.However the intrinsic influence of grain shapes on R might be more difficult to understand given that natural variations in shape are generally coupled with changes in grain size.No general trend is obvious with our limited set of samples and this point should be studied further.
The double peak that appears for R in the 1.5 µm absorption band for S1 (dendritic crystals) in Fig. 6a, also appears in the SnowRAT simulations for randomly oriented cylinders (Fig. 9).Both observations indicate that the double peak may be caused by elongated forms or faceted crystals (dendritic Atmos.Chem.Phys., 10, 2507-2520, 2010 www.atmos-chem-phys.net/10/2507/2010/crystals, cylinders, columns ...).Simulations have been done to investigate the origin of this double peak using cylinders.
Even with artificial increase of the imaginary refraction index of ice, the two maxima in the BRDF (at limb and at (θ v =30 • , φ=180 • )) remain.This indicates that the double peak may solely involve single or multiple reflections at the surface of the grains without transmission through the grains.
The double peak has rarely been measured in the past since angular field of view of the sensors is typically too large.In addition, the incident zenith angles are generally higher than 40 • .

Model/measurements discrepancies
The rainbow that appears in Fig. 8 is due to the fact that perfect spheres are used for this simulation.The phase function of spheres (Fig. 2, Mishchenko et al., 1999) presents a local maximum at 138 • scattering angle due to internal reflections.This explains the rainbow when the sun illuminates droplets.The higher the absorption is, the more contrasted the rainbow is; this point is confirmed in Fig. 8.This crescent or rainbow -which is not observed for natural snowpacks, is absent for any other grain shape used in the simulations.
Some other discrepancies still remain between the models and the measurements.The darkening at grazing angles is stronger in the model results than in our measurements.It is partly due to the fact that forward scattering is usually overestimated using spheres in models especially at near infrared wavelengths (Jin et al., 2008).At wavelengths shorter than 1 µm, the discrepancies are mainly explained by the fact that a real snowpack (sample) has surface roughness whereas models assume a smooth surface (Jin et al., 2008).At wavelengths longer than 1 µm, a possible explanation of the slight differences between measured and modelled R is that our samples are a complex mix of shapes and sizes while the model considers only one single shape with power law size distribution.

Conclusions
This paper presents a large set of direct measurements of BRDF for different types of snow.The comparison with modelled BRDF and results in literature allow to explain the main BRDF variations as a function of viewing and lighting angles, wavelength, size and shape of grain.
The first point to underline is that the variations of the anisotropy factor with wavelength are controlled by the ice absorption coefficient.For wavelengths shorter than 1 µm, the most noteworthy effect at near vertical incidence is the darkening at grazing angles.This effect is a consequence of dominant multiple scattering within the snowpack.In contrast, for wavelengths longer than 1 µm and/or large incident zenith angles, forward scattering is stronger because absorption is high and single scattering prevails and thus the anisotropy factor is mostly controlled by the phase function.Grain size and shape have a great influence.However, their respective effect on the anisotropy factor are difficult to predict at near IR wavelengths.For a given shape, the anisotropy increases with increasing size but comparison between snow grains with different shapes and sizes are more complicated.For elongated or faceted shapes such as dendritic crystals, columns or cylinders, two maxima appear on the anisotropy factor patterns.
Photon tracing and radiative transfer models predict anisotropy factors in general close to measurements.Using non-spherical shapes allows to simulate feature as the double peak, avoid artefact such as rainbow that appears for spheres and probably contribute to better agreements between models and measurements at 30 • incident zenith angle for 1.0 and 1.03 µm (Fig. 8, supplement: http://www.atmos-chem-phys.net/10/2507/2010/acp-10-2507-2010-supplement.pdf).Some discrepancies still exists most likely due to the complex mix of different crystal shapes and to the surface roughness of natural snowpack.
These results allow to estimate the error implied while considering snow as a Lambertian surface for processing remote sensing data.Furthermore, these results make accute retrieval of snow surface spectral albedo from remote sensing reflectance data possible.

Fig. 3 .
Fig. 3. S3 Anisotropy factor, R(λ) for different viewing angles in the principal plane.Incident angle is 30 • .Vertical lines are the same as in Fig. 2. Gray areas correspond to wavelengths where photometric accuracy is reduced due to absorption by H 2 O vapor or inaccurate calibration of the reference.Two polar plots in the principal plane, with all the measured viewing angles have been added at 0.7 and 1.7 µm above the chart.

AtmosFig. 4 .
Fig.4.S3 anisotropy factor, R(θ v ,φ), at 0.6 and 1.5 µm for zenith incident angles 0 • , 30 • and 60 • .The polar angle corresponds to the relative azimuth, φ, between the viewing and the incident azimuth and the polar radius to the viewing angle, θ v .The incident beam comes from the right and the forward direction is toward left.The three circles inside each plot represent viewing zenith angles of 30 • , 60 • and 70 • .The crosses show the measurements used to generate the isolines.The spectral albedo used to calculate R is indicated below each chart.

Fig. 5 .
Fig. 5. Ratio of the anisotropy factor R between samples S1 and S3 R(S1) R(S3) at 0, 30 and 60 • incident beam and for two wavelengths.The spectral albedo used to calculate R for S1 is indicated below each chart.S3 spectral albedo is indicated in Fig. 4.

Table 1 .
Summary of physical properties of snow samples.
a Initial means measurement as the sample is collected i.e. before storage in cold rooms.