Analysis of snow bidirectional reflectance from ARCTAS spring-2008 campaign

A. Lyapustin, C. K. Gatebe, R. Kahn, R. Brandt, J. Redemann, P. Russell, M. D. King, C. A. Pedersen, S. Gerland, R. Poudyal, A. Marshak, Y. Wang, C. Schaaf, D. Hall, and A. Kokhanovsky University of Maryland Baltimore County, Baltimore, MD, USA NASA Goddard Space Flight Center, Greenbelt, MD, USA University of Washington, Seattle, Washington, USA Bay Area Environmental Research Institute (BAERI), Sonoma, CA USA NASA Ames Research Center, Moffett Field, CA, USA University of Colorado, Boulder, CO, USA Norwegian Polar Institute, 9296 Tromsø, Norway Science Systems and Applications, Inc., Lanham, MD, USA Boston University, Boston, USA Institute of Environmental Physics, University of Bremen, 28359 Bremen, Germany

for the glint region (azimuthal angles ϕ < 40 • ), the best fit MRPV and RTLS models fit snow BRF to within ±0.05.The plane-parallel radiative transfer (PPRT) solution was also analyzed with the models of spheres, spheroids, randomly oriented fractal crystals, and with a synthetic phase function.The latter merged the model of spheroids for the forward scattering angles with the fractal model in the backscattering direction.The PPRT solution with synthetic phase function provided the best fit to measured BRF in the full range of angles.Regardless of the snow grain shape, the PPRT model significantly over-/underestimated snow BRF in the glint/backscattering regions, respectively, which agrees with other studies.To improve agreement with experiment, we introduced a model of macroscopic snow surface roughness by averaging the PPRT solution over the slope distribution function and by adding a simple model of shadows.With macroscopic roughness described by two parameters, the AART model achieved an accuracy of about ±0.05 with a possible bias of ±0.03 in the spectral range 0.4-2.2µm.This high accuracy holds at view zenith angles below 55-60 • covering the practically important range for remote sensing applications, and includes both glint and backscattering directions.

Introduction
Due to its high reflectance, snow is one of the key factors defining the surface radiative budget in polar regions and affecting global climate of the Earth.The albedo of thick snowpack is governed primarily by the snow grain size and level of impurities, such as soot or dust, deposited from the atmosphere.Monitoring snow properties from space requires an accurate knowledge of its bidirectional reflectance factor (BRF) (Stroeve and Nolin, 2002;Stroeve et al., 2005;Stamnes et al., 2007;Scambos et al., 2007;Painter et al., 2009;Lyapustin et al., 2009).Accurate modeling of snow reflectance is equally important for atmospheric remote sensing of cloud properties and vital for aerosol retrievals over snow-covered regions, which presently remains a largely unresolved problem.
Beginning with the classic work of Wiscombe and Warren (1980), snow reflectance has been studied extensively using a plane-parallel radiative transfer model.The accumulated body of measurements and modeling efforts converges in understanding that the bidirectional reflectance of snow is less anisotropic than predicted by the RT model (e.g., Warren et al., 1998;Painter and Dozier, 2004;Hudson et al., 2006;Hudson and Warren, 2007).The main errors appear around the principal plane, where the model significantly overestimates reflectance for forward scattering angles and underestimates it in backscattering directions.Because of the apparent nonsphericity of snow grains, a number of studies investigated the effect of snow particle shape on the modeled reflectance (Mishchenko et al., 1999;Xie et al., 2006;Jin et al., 2008).Yang and Liou (1998), Kokhanovsky and Zege (2004), and Jin et al. (2008) found that using the nonspherical model of snow grains improves agreement with measurements.The latter work, which included variable particle shapes and particle surface roughness, compared modeled results with tower measurements from 32 m made over Antarctica (Hudson et al., 2006).Although microscopic roughness (with scale much less than the grain size) smoothed the difference, the best theoretical predictions still significantly overestimate tower snow reflectance in the forward scattering angles and underestimate it in the backscattering region.Overall, the authors achieved agreement between modeled radiance and tower measurements within ±10% for view zenith angles ≤60 • , with the abovementioned asymmetry between the forward and backscattering directions.These studies help clarify the point that the microphysical variability alone cannot fully capture all features of snow reflectance, and macroscopic (∼10 cm) effects of non-flat snow surfaces and shadows need to be incorporated into the model to achieve closure with measured snow BRF, even at the tower-scale footprint, not to mention the satellite footprint.Hudson and Warren (2007) have demonstrated the effect of snow surface roughness by contrasting snow reflectance under clear skies with that from clouds over a snow at different view zenith angles.Nolin et al. (2002) used multi-angle MISR data to characterize roughness of sea ice and ice sheets over Greenland and Antarctica.A study by Leroux and Fily (1998) showed that modeling the effect of sastrugi improved agreement with snow BRF measurements over Antarctica.
In this work, we study the effect of macroscopic surface roughness, including the slope distribution of reflecting facets and shadows cast by sastrugi, using a simple model.We also study the accuracy of the common analytical BRF models used in operational satellite data processing, including the reciprocal Ross Thick -Li Sparse (RTLS, Lucht et al., 2000) and Modified Rahman-Pinty-Verstraete (MRPV, Martonchik et al., 1998) models.The three-parameter RTLS and MRPV models are used in the MODIS and MISR processing, respectively.A specialized Asymptotic Analytical Radiative Transfer (AART, Kokhanovsky and Zege, 2004) model, which has been actively explored for the snow grain size retrievals (Zege et al., 2008;Tedesco and Kokhanovsky, 2007;Lyapustin et al., 2009), is also evaluated in this paper.In parallel, we assess the atmospheric correction error over snow due to the Lambertian assumption, which is the basis of MODIS operational land processing (Vermote and Kotchenova, 2008).
Our analysis uses airborne measurements obtained by the Cloud Absorption Radiometer (CAR, King et al., 1986;Gatebe et al., 2003) during the ARCTAS Spring campaign (Jacob et al., 2009).The atmospheric correction algorithm uses ancillary aerosol measurements from the Aerosol Robotic Network (AERONET, Holben et al., 1998) and data from the Ames Airborne Tracking Sunphotometer (AATS, Russell et al., 2005Russell et al., , 2007) ) to derive snow BRF.CAR provides comprehensive spectral coverage from the UV through shortwave infrared spectral region, and has unprecedented angular resolution for an airborne sensor of 1 • in zenith and azimuthal angles.Several representative experiments were conducted at Barrow, Alaska, over Elson Lagoon.Different flight altitudes allow us to test the atmospheric correction algorithm and get insight into snow spatial homogeneity at scales of ∼0.2-2 km.The atmospheric correction algorithm is independently evaluated against surface measurements of snow albedo.
This paper is organized as follows.Section 2 provides a brief description of the NASA P-3B aircraft flights, focusing on CAR measurements, and Sect. 3 describes ground albedo experiment.The atmospheric correction algorithm is described in Sect.4, along with its evaluation.Based on derived snow BRF, Sect. 5 examines analytical BRF models.Finally, macroscopic surface roughness is explored in Sect.6, followed by Conclusions.The analytical BRF models and the model of surface roughness are described in the Appendices.

Description of experiments
The ARCTAS Spring campaign was conducted in 1-21 April 2008, with the aim of studying physical and chemical processes in the Arctic atmosphere and related surface phenomena, as part of the International Polar Year.NASA deployed two aircraft from Fairbanks, Alaska, the DC-8, instrumented primarily for atmospheric chemical sampling, and the P-3B, which carried a payload, including the CAR and AATS instruments, designed to study aerosols and the radiation environment (Jacob et al., 2009).A third NASA aircraft, the B200, was stationed at Barrow, Alaska, and carried the NASA Langley High-Spectral-Resolution Lidar (HSRL;Hair et al., 2008).In addition, a number of AERONET and AEROCan sun and sky-scanning photometers (Holben et al. 1998;Bokoye et al., 2001) operated in the study region during this period, including one at Barrow and one at Eureka, Canada.
The ARCTAS directional snow reflectance experiment addressed multiple objectives, as the albedo and snow reflectance angular variation are critical for surface characterization, energy balance calculations, and the remote sensing of atmospheric properties, including aerosol amount and type.For this experiment, we collected coincident observations of the snow-covered surface from platforms at multiple levels.The suite of measurements included ground-based Analytic Spectral Device (ASD) FieldSpec radiometer observations and direct snow samples, airborne instrument data, including CAR and AATS remote sensing from the P-3B, and near-coincident observations from the NASA Earth Observing System's Terra satellite Multi-angle Imaging SpectroRadiometer (MISR) and MODerate-resolution Imaging Spectroradiometer (MODIS).The primary target site for this experiment was Elson Lagoon (71.3 • N, 156.4 • W), near Barrow, Alaska, which was studied by the P-3B payload on 7 and 15 April.The HSRL joined the P-3B on 7 April, and coincident surface measurement took place the week of 15 April.Other sets of measurements over snow surface was made with the P-3B instruments on 8 April near Eureka, Canada (80.5 • N, 90.2 • W) and on 9 April near Axel Heiberg Island, Canada (79.9 • N, 100.9 • W).The present study focuses on the Elson Lagoon events.
3 Albedo measurements at the surface A site was selected 10 km east (upwind) of Barrow, Alaska, on Elson Lagoon.This lagoon is a protected arm of the Beaufort Sea.The surface consisted of flat land-fast first-year sea ice of thickness ∼1.5 m, covered with 25-40 cm of snow with density 0.35 g cm −3 .The snow surface roughness was dominated by small sastrugi, with typical height ∼5 cm and spacing ∼1 m (Fig. 1).The lagoon is a homogeneous area large enough to ensure a uniform representative footprint in the field of view of the CAR from the P3 airplane.Spectral albedo was measured using a FieldSpec Pro JR spectroradiometer manufactured by Analytical Spectral Devices, Inc. (ASD) (Kindel et al., 2001).The instrument records radiance every 1 nm from 350 to 2500 nm wavelength, with 3-to 30-nm spectral resolution.Two ASD instruments were deployed at Barrow, one from the University of Washington (UW) and the other from the Norwegian Polar Institute (NPI).A "cosine collector", receiving radiation from a hemisphere, is supported by a 1.6-m rod, mounted on a tripod.The rod has bubble levels at both ends, so the cosine collector can be rotated to view the upward and downward hemispheres alternately.The light received by the cosine collector is directed by a fiber-optics guide to the ASD detector.Although the two radiometers were identical, their cosine collectors were not.The UW cosine collector employs an isotropic spectralon reflecting plate while the NPI instrument measures transmitted light through a diffuser type cosine collector.The cosine collectors were calibrated as in Grenfell et al. (1994).The UW instrument had better cosineresponse than the NPI instrument at infrared wavelengths, and the reverse was true at visible wavelengths; this affects our selection of data below.
Under overcast skies with diffuse incidence, the cosine collector, support rod, tripod, radiometer box, computer, and sled all intercept some of the downward solar radiation.To account for this reduction, a "shadowing correction" based on geometric analysis was applied to the albedo measurements.For diffuse incidence the shadowing correction was estimated as 1.7% for the UW instrument and as 0.6% for the NPI instrument.For direct incidence with low sun the shadowing correction is smaller.
Spectral albedo is usually measured more accurately under overcast sky when both upward and downward radiation is diffuse and errors caused by leveling the instrument, by deviations from cosine response, and by non-horizontality of the surface are much smaller than under the direct-beam incidence.The surface spectral albedo was measured on six days in mid-April 2008 under both clear and overcast skies including 15 April when simultaneous measurements were made with the CAR from the NASA P-3B aircraft.The surface air temperature during the surface albedo experiments was in the range −14 to −21 • C throughout 14-17 April.On the 19th it was warmer, −1.3 • C.

Measurements on 15 April
Fresh snow had fallen on 14 April with minimal redistribution by the wind.Albedo was measured on 15 April under clear skies conditions at five sites near the center of a designated clean area along a westerly track at 50-m intervals.Variation in measured albedo at a given site was significant due to slight leveling errors in positioning the cosine collector.These errors were corrected using surface albedo measured five days later under overcast sky where the diffuse conditions minimize errors due to leveling and cosine response.The leveling-error correction was applied as a spectrally-constant scale factor, chosen so as to adjust the 490-nm albedo to equal that measured on 19 April under overcast cloud.The measurements on 19 April were not affected by leveling errors because the incident radiation was diffuse, and the high albedo at 490 nm is insensitive to grain size variation (Wiscombe and Warren, 1980, Fig. 8a).Furthermore, the solar zenith angle during the 15 April measurements was ∼60 • , which is close to the effective incidence angle of diffuse irradiance, so the effect of zenith-angle difference on albedo of the two days should be negligible, particularly at 490 nm where albedo is insensitive to zenith angle (Wiscombe and Warren, 1980, Fig. 11).The leveling-error correction at 490 nm was less than 1.5%, while the agreement between the NPI and UW at this wavelength was better than 1%.
A correction was applied for the radiometer and support shadows, which blocked a small portion of the downward looking field-of-view.Due to differences in cosine collector design, the NPI ASD was more accurate for wavelength below ∼1 µm, whereas the UW ASD provided a better signal/noise ratio beyond ∼1 µm.For this reason the UW albedos were scaled by a factor of 0.986 at all wavelengths to match the NPI albedo in the 800-850 nm band where both instruments performed well.The final result (Fig. 2) is a combination of the NPI albedo at UV and visible wavelengths and the UW albedo in the near-infrared.
After corrections and scaling, variation among the five sites was significantly greater than the variation of repeated measurements at a particular site.The measured albedos at each site were averaged; the final albedo represents an average among 5 sites.The error envelope in Fig. 2 includes variation between sites, estimated error in the shadowing correction and error in assigning the 490-nm albedo.4 Processing CAR data

Atmospheric correction algorithm
To process CAR data over snow, a specialized atmospheric correction algorithm has been developed.It is based on an accurate semianalytical expression for the atmospheric radiance derived with the Green's function method (Lyapustin and Knyazikhin, 2001;Lyapustin and Wang, 2005).When used with the RTLS model (see Appendix A), this solution makes it possible to express upward reflectance at an arbitrary altitude z in the atmosphere as an explicit function of RTLS model parameters where R D is the atmospheric (path) reflectance, F k (k= L, G, V ) are functions of view geometry and atmospheric properties, µ 0 and µ are cosines of the solar and view zenith angles, respectively, and ϕ is the relative azimuth angle.Functions F k have a weak dependence on surface reflectance through a multiple reflection factor, α = (1 − q(µ 0 )s) −1 , where q is surface albedo and s is spherical albedo of the atmosphere.R nl is a small term that is nonlinear in surface reflectance.
The quasi-linear form of Eq. ( 1) leads to a very efficient iterative minimization algorithm for the root-mean-square error for three parameters of RTLS model K: In this expression, R is a measured CAR reflectance, the index j denotes different CAR view angles, and n is the iteration number.Although it is not shown explicitly, functions F k are also updated every iteration, according to the current value of parameter α.Equation ( 2) provides an explicit leastsquares solution for coefficients K: where In the first iteration, the non-linear term is computed for an assumed spectrally dependent albedo, for example q (0) (vis) = 0.6, or q (0) (2.2 µm)=0.05.Convergence is achieved in 4-5 iterations over bright snow in the CAR blue band (0.48 µm), and it takes fewer iterations at longer wavelengths where snow is less reflective.
Once the RTLS coefficients are computed, the diffuse component of the surface-reflected signal (R Dif ) is calculated at the altitude z of the P-3B flight, and the snow bidirectional reflectance factor ρ, representing specific measurements, is computed from the direct reflected term.For this purpose, we single-out the direct term and re-write Eq. ( 1) as follows: Here, τ 0 is the column optical depth through the atmosphere at the wavelength of interest.Finally, once the snow BRF is computed, the best-fit parameters of the MRPV and AART models are retrieved and the rmse is evaluated.The MRPV parameters are computed using logarithm transformation, following the MISR surface reflectance algorithm (Martonchik et al., 1998).The required RT calculations are performed with the SHARM code (Lyapustin and Wang, 2005).Assuming all the ozone is located above the aircraft, the CAR data are preliminarily corrected for ozone absorption along the solar path.

Ancillary parameters
Table 1 lists the CAR spectral bands, the in-band solar irradiance used to convert digital numbers into reflectance units, and the assumed optical thickness of column gaseous absorption.The latter includes 0.4 cm column water vapor that is close to the AATS above airplane value at the lowest flight altitude (Schmid et al., 2003;Livingston et al., 2008) and to the AERONET column measurement at the Barrow site during these experiments.The band-average gaseous absorption, τ g , is computed from monochromatic values, τ g (λ), using the equation: where h λ is the CAR relative spectral response function, and F λ is the exoatmospheric spectral solar irradiance.The optical thickness of absorbing gases is calculated for a carbon dioxide concentration of 340 ppm, and other gas concentrations corresponding to the sub-Arctic Winter atmospheric model (Kneizys et al., 1996).The calculations included absorption of six major atmospheric gases (H 2 O, CO 2 , CH 4 , NO 2 , CO, N 2 O) with line parameters from the HITRAN-2000 (Rothman et al., 2003) database, using the Voigt spectral line shape and the Atmospheric Environmental Research continuum absorption model (Clough et al., 2005).Because of the water vapor, oxygen and other gas absorption band structure, the effective in-band transmittance of the atmospheric column changes with air-mass approximately as √ m rather than ∼ m.For this reason, Table 1 shows two numbers for effective column optical thickness, one for vertical path, and another for a solar zenith angle (θ 0 ) of 67.2 • .The second value, which was used in the processing algorithm, gives a correct transmittance for the incident radiation for a given θ 0 , although it may introduce some view-angle dependent error on the light path from the surface to the aircraft.This error, however, is small because the flights were at low altitude (0.2-1.7 km) and the gaseous absorption is low.Table 1 also lists the in-band optical thickness of ozone absorption for 445 DU that was observed over Barrow during 7-15 April by the Ozone Monitoring Instrument (OMI).To verify CAR calibration in reflectance units, we provide irradiance F λ obtained by integrating data from Solar Irradiance Monitor on SORCE (Harder et al., 2000) collected on 8 April.
The aerosol optical thickness above the P-3B aircraft was measured by the AATS onboard the aircraft, and the total column aerosol optical thickness was acquired by the AERONET sunphotometer at Barrow.The column water vapor provided by AERONET was low, ∼0.4 cm.Table 2 summarizes the P-3B flights with CAR measurements over snow, including dates, average height, solar zenith angle, and AERONET and AATS aerosol optical thickness.

Algorithm and CAR data analysis
Snow is the most reflective of the common land surface types, having an albedo close to unity in the visible spectrum, and data analysis over pure snow raises accuracy standards for the processing model as well as for sensor calibration accuracy.It may be mentioned in advance that CAR has provided excellent spectral BRF data for snow anisotropy analysis.However, we encountered a problem regarding the overall reflectance magnitude at the red band (0.68 µm) in the experiment conducted on 7 April, which is described and analyzed in this section.
The atmospheric correction algorithm developed here automatically computes surface albedo once the RTLS model parameters are found.The RTLS model approximates snow BRF rather well, with standard deviation in the blue band below 3% relative to nadir.Forward computations with derived snow BRF reproduce the CAR measurements with an accuracy of better than ±0.5%.In addition, both spectral BRF and albedo derived from CAR flights at different altitudes on 7 April agree with each other to better than ∼0.3%, and also agree with results obtained over the same location on 15 April at lower solar zenith angle.These analyses support the conclusion that the algorithm works correctly.However, the albedo derived for the CAR red band on 7 April was found to be slightly but systematically higher than unity, which violates energy conservation.The derived spectral albedo was compared with surface albedo measurements collected on 15 April.Figure 2 shows the envelope of measured albedo from minimal to maximal values over the five sites; they are within about ∼2% in the visible bands.The triangles and circles show albedo derived at the CAR central wavelengths on 7 and 15 April, respectively.The overall agreement between CAR and ground albedo is good, especially given possible snow inhomogeneity, and the difference between the CAR footprint and the coverage of the ground measurements.Except in the red (0.68 µm) and NIR (1.27 µm) bands, the retrievals are within or very close to the envelope of measured albedo.In the red band, however, derived albedo is about 4-5% higher than the maximal measured value.It exceeds unity on 7 April (q = 1.015) and is very close to unity (0.999) on 15 April.
We analyzed the main factors affecting the CAR red band results.
1) To exclude the possibility of error in the calibration conversion coefficients from radiance digital number (DN) counts to reflectance, the irradiance F λ used in the CAR calibration was compared with irradiance from the Solar Irradiance Monitor (Harder et al., 2000) integrated over the CAR spectral response (see Table 1).Although the use of SIM irradiance reduces reflectance in the red band, the effect is only about 0.1%.
2) Because of the relatively large solar zenith angle, it is important to accurately model absorption by well-mixed gases.The CAR red channel overlaps with the weak oxygen B-band, the main absorber in this channel apart from ozone.Our current modeling of absorption relies on the HITRAN-2000 edition.Although the new HITRAN-2008 database has been released, there has been only a minor change in oxygen line parameters in the region of interest (Rothman et al., 2009).It cannot explain our result because the absorption optical thickness would need to decrease by approximately 0.015 (or by ∼80%) to reduce the albedo by 4-5%.
3) Two algorithm-related factors were investigated as part of this analysis.First, although CAR provides the full hemisphere of upward directions, the maximal view zenith angle (θ ) in our processing is limited to about 75 • to conform to the limits of the plane-parallel radiative transfer model.To assess whether the upper θ limit impacts derived RTLS coefficients and albedo, θ max was varied from 60 • to 85 • .The solution was very stable, with maximal change below 0.2% in the blue band and negligible in the other bands.Second, the accuracy of the best-fit RTLS model is not uniform on the hemisphere of upward directions, and potentially, high outliers may bias the derived albedo.To investigate this factor, we added a separate albedo retrieval algorithm that does not depend on the model of surface BRF.It consists of three steps: in the first step, CAR measured radiance is directly integrated using Simpson's quadrature to obtain upward flux at the flight altitude z, , where F up Th is a "theoretical" flux computed with the best-fit RTLS parameters and given atmospheric profile.Although individual "theoretical" fluxes at levels z and z = 0 depend on parameters of the BRF model, this dependence disappears in the ratio.Finally, surface albedo is computed by dividing obtained reflected flux F up (0;µ 0 ) by the "theoretical" incident flux F dn Th (0;µ 0 ).We found that this approach generally decreases the albedo by 0.8% in the blue band, and by 0.5-0.2% at longer wavelengths, the difference diminishing with wavelength.This albedo reduction is taken into account in the data shown in Fig. 2.This analysis shows that the observed differences in the red band cannot presently be explained by known algorithmrelated factors, and one plausible explanation is a calibration uncertainty of CAR.The radiometer was calibrated using a calibration sphere, which may have an absolute uncertainty within 5%.Given this source of uncertainty, the bias in the red band has yet to be explained, particularly in view of the absence of such a bias in the nearby blue band.Our analysis will continue with further investigation of issues related to the processing algorithm and absolute CAR calibration.

Bidirectional reflectance of snow
Figure 3 shows the derived snow BRF in the CAR red band.The first three images show results for 7 April from three different observing altitudes, 1.7, 0.64, and 0.2 km, and the last image shows retrievals for 15 April for measurements from 184 m altitude.The BRF has a very consistent pattern, with a shape dominated by two features -the increase of reflectance in the glint (forward scattering) direction, and a smaller reflectance increase in the backscattering direction.The dark feature in the hotspot direction in the last two images (c, d) corresponds to the airplane shadow.
The radiometer scans different surface points while obtaining the total BRF.The CAR footprint decreases, and spatial resolution increases, at lower flight altitudes; the effect of surface inhomogeneity becomes clearly visible at altitudes of ∼0.2 km and below, on both 7 and 15 April.Analysis of BRF at angles close to nadir, for example at θ 3 • , makes it possible to evaluate snow reflectance variability (inhomogeneity) as ∼0.05 for flight altitudes above 0.6 km (>10 m footprint at nadir) and ∼0.15-0.2 for 0.2 km altitude (4 m nadir footprint).
Snow BRFs remain close in the blue-NIR range, 0.47-0.87µm.Although the imaginary part of the ice refractive index increases with wavelength, it remains very low in this spectral range, as does snow absorption.In these conditions, BRF is dominated by multiple scattering, and BRF anisotropy is low.At longer wavelengths, absorption becomes progressively more prominent, and as the magnitude of snow reflectance drops, reflectance anisotropy increases.This result is shown in Fig. 4, which shows the spectral dependence of snow BRF in the principal plane and in the cross-plane for flight segment 2001b (0.64 km).The blue band saturates around the maximum glint region (θ > 67-71 • ) which explains decrease of the measured CAR signal and the retrieved BRF.Taking the ratio Ratio =ρ(θ=60 • )/ρ(θ=0 • ) at ϕ=0 • in the principal plane as a measure of anisotropy, one can see the increase of anisotropy with wavelength quantitatively: Ratio={1.44;1.48; 1.69; 1.93; 1.99; 3.66; 3.7} at λ={0.68; 0.87; 1.03; 1.22; 1.27; 1.65, 2.2}.Here, we omitted the blue band because of the saturation problem.
The cross-plane (ϕ=90 • ) reflectance increases systematically with the view zenith angle.Although anisotropy of reflectance is notably lower than in the principal plane, nonetheless snow reflectance is significantly non-Lambertian.

Analysis of MRPV and RTLS models
The RTLS and MRPV models (see Appendix A) are used in the MODIS (Schaaf et al., 2002) and MISR (Martonchik et al., 1998) land algorithms, respectively, including snowcovered surfaces.The accuracies of these models were studied experimentally over different land cover types (Privette et al., 1997;Bicheron and Leroy, 2000), however we are not aware of such analysis over snow.The CAR dataset presents an excellent opportunity to perform an accuracy analysis over pure snow.In order to remedy a possible bias in the CAR red band discussed in Sect.4.3 and shown in Fig. 2, the derived BRF was renormalized by matching retrieved red band albedo with the mean value from ground-based measurements.This renormalization does not affect the reflectance angular distribution.The retrieved and normalized snow BRF is used below to evaluate goodness of fit by the RTLS and MRPV models.
The results of model analysis are shown in Fig. 5 for the red (0.67 µm) and NIR (1.22 µm) bands that are often used in snow property retrievals.The left image shows the retrieved BRF followed by images showing accuracy of the RTLS and MRPV models.The accuracy, computed as (ρ − ρ Model ), is shown for the full range of values and for a reduced range of values (±0.05).One can see that though both models show a relatively similar performance, differences are clearly visible near the principal plane.The best-fit MRPV model provides a better approximation in the backscattering directions near the principal plane, where the RTLS model significantly underestimates reflectance.On the other hand, the best fit RTLS model approximates the reflectance at forward scattering angles much better: the maximal error of the RTLS model is 0.12-0.18,whereas the MRPV model error reaches 0.28-0.45 in the glint direction.
Figure 6 shows the error of atmospheric correction at wavelengths 0.47, 0.67, 0.87 and 1.22 µm when a Lambertian surface model is assumed.As is well known, the Lambertian assumption reduces the anisotropy of the BRF.It overestimates BRF at low view zenith angles, where BRF is lowest, and significantly underestimates it where the BRF is high, in the glint and hotspot regions.The error decreases with increasing wavelength.For example, the average reflectance overestimation at low θ is 0.02-0.025,0.01-0.015,0.008-0.012,and 0.006-0.008at 0.47, 0.67, 0.87, and 1.22 µm, respectively.The maximum error at θ=60 • in the forward scattering direction at the same wavelengths is −0.18, −0.10, −0.064 and −0.026, respectively.In all the ARCTAS snow experiments, the atmosphere was relatively clear (see Table 2).The error will be higher under hazy conditions, due to increased atmospheric scattering.

Analysis of AART model
The Analytical Approximate Radiative Transfer model (see Appendix A) predicts snow bidirectional reflectance as a function of snow grain size and level of absorbing impurities, such as soot.The analytical form of this model makes it very convenient for snow grain size retrievals, and several approaches have recently been developed (Tedesco and Kokhanovsky, 2007;Lyapustin et al., 2009;Kokhanovsky and Schreier, 2009).The AART model was tested against snow pit measurements ( Kokhanovsky et al., 2005).Contrary to our results, the AART model was found to overestimate measured BRF in the glint region under specific conditions of measurements.Nevertheless, this work, and a broadscale analysis of MODIS data over Greenland (Lyapustin et al., 2009), showed that the model is not valid in forward scattering directions near the principal plane (at ϕ < 40 • ).Kokhanovsky et al. (2005) found that it may produce large errors at angles close to the backscattering direction as well.These errors, discussed in the introduction, are inherited from the plane-parallel radiative transfer solution (function R 0 in Eq. (A3) of Appendix A).
In this section, we study the effect of macroscopic surface roughness on snow bidirectional reflectance.Roughness affects snow reflectance in two ways, via the variety of slopes, which requires averaging the plane-parallel solution over a slope distribution function, and via shadows.In this work, we use a Gaussian slope distribution function with averaging procedure described in Appendix B. The model of snow shadows described in Appendix C may be considered as oversimplified, and yet it allowed us to achieve a relatively good match with BRF derived from the CAR measurements.The resultant model of snow surface roughness depends on two parameters, standard deviation of slope distribution (σ ) and sastrugi density (D).

Averaging over slope distribution function
Figure 7a-f shows the effect of averaging over a slope distribution function for several values of standard deviation σ =0, 0.2, 0.4.Figure 7a-c shows results for a model using spheroids (Dubovik et al., 2006) to compute the snow scattering function.The reflectance of a semi-infinite medium is computed with the SHARM code (Lyapustin, 2005) for zero absorption and an optical thickness of 20 000, which provides an asymptotic reflectance limit.The results are shown for a grain size of 50 µm.Because there is no absorption, we did not find any noticeable difference in the radiative transfer solution when the grain size varied from 20 to 200 µm.Figure 7a shows the deviation of the planeparallel solution (σ =0) from the experimental BRF computed as R meas − R model .It overestimates the CAR BRF by 1.21 in the glint region maximum, and underestimates it by 0.31 in the backscattering direction.Averaging over slopes decreases the forward peak of the theoretical reflectance and increases its backscattering value, thereby reducing the differences with experimental data.At σ =0.4, the difference in the forward and backscattering peaks is reduced by a factor of 3 (−0.42)and 2 (0.166), respectively, at θ=76 • .
Figure 7d, e shows comparisons for an averaged BRF (σ =0.4) computed with a Mie model (50 µm spheres) and a model of randomly oriented fractal crystals (Mishchenko et al., 1999).In addition to errors in the rainbow angles, which are usually not observed over snow, at least for visible wavelengths, the Mie model fits the CAR BRF poorly in the forward scattering region, where the error remains inhomogeneous and large over a wide range of azimuths and view zenith angles.The fractal model has an interesting pattern around the forward scattering peak, underestimating BRF at its center and overestimating it at angles 30-60 • off center.For the experiment geometry (θ 0 =68.6 • ), this suggests that the peak of the fractal phase function is too narrow and the phase function is underestimated at scattering angles below ∼50 • , while at the same time, it is overestimated at scattering angles ∼60-90 • .Comparing phase functions for the spheroid and random fractal models, Fig. 8 shows that this is indeed the case.Except for the extremely narrow diffraction peak, the fractal phase function is significantly lower at angles below 28 • , and it is higher than the spheroidal model at angles above ∼39 • .Overall, the fractal model gives a much better fit in the backscatter region, where the error is a factor of ∼2 lower than that of the spheroid model.Given these considerations, we constructed a synthetic phase function from the 50 µm spheroid (SPD) and fractal models as follows: x Synth (γ )=x SPD (γ ) at scattering angles γ < 37 x Synth (γ ) = x Fract (γ ) at γ ≥ 90 • ; and finally, in the range 60 • < γ < 90 • , x Synth (γ ) was obtained by linear interpolation between the x Synth (60 • ) and x Synth (90 • ) values.For normalization, the resulting scattering function was reduced by a factor of 1.103.Figure 7f shows this result for the synthetic BRF with slope averaging (σ = 0.4).One can see the reduction in the total difference with the CAR BRF down to (−0.1; 0.12) for the full range of θ ≤ 76 • , which is an improvement over all other cases.It is possible that specialized particle shapes modeling ice crystals (e.g., Xie et al., 2006;Jin et al., 2008) may provide an even better agreement with experimental data.

Effect of shadows
In contrast to the slope averaging, shadows are not visible in the backscatter direction, where they do not change the absolute value of reflectance.In our case, this represents the angles θ ≥ θ 0 at ϕ = 180 • .On the other hand, shadows act to increase the relative value of the backscattering peak by decreasing reflectance at other angles.Figure 7g-l shows the effect of shadows on the BRF for two values of increased sastrugi density (D = 0.04, 0.08).The images (g-i) are computed for the spheroid model at σ = 0.4, and the images (j-l) show results for the synthetic model at σ = 0.3.Figure 7i and l reproduces Fig. 7h and k at an expanded scale to show the fine details.Here, the black colors show values above the upper limit and below the lower limit of the color bar.One can see that the best approximation of the CAR BRF for flight 2001b is obtained with the Synthetic model, providing accuracies within −0.02 to +0.07 at θ ≤ 60 • .Finally, Fig. 7m shows the theoretical model fit for experiment 2005a, which had a different solar zenith angle (61.95 • ).In this case, a good approximation, within ±0.06 at θ ≤ 60 • , was achieved with the Synthetic model and parameters σ = 0.2, D = 0.10, close to values for the previously considered experiment 2001b.

Spectral analysis
If snow reflectance R 0 is known for the nonabsorbing shortwave channel, then the AART model can predict the spectral dependence of snow BRF analytically (Eq.A3, Appendix A): Assuming the function R 0 is computed with the synthetic aerosol model and parameters (σ = 0.3, D = 0.08), we analyzed the accuracy of the AART model at 0.681, 0.871, 1.03, 1.22, 1.27, 1.654 and 2.204 µm wavelengths, for CAR flight 2001b.Based on the theoretical derivation, the AART model is valid for 4π χ d λ 1, limiting it to the visible-NIR spectral range of low absorption.So, our analysis at 1.654 and 2.204 µm wavelengths for high snow absorption should be considered as formal only.The imaginary part of the refractive index of ice was taken from Warren and Brandt (2008).
Conducting this analysis requires knowledge of parameter d.The snow grain size (diameter) for this exercise was first evaluated from CAR BRF data using a method developed for MODIS (Lyapustin et al., 2009), and based on the ratio of the NIR (λ 1 = 1.22 µm) to red band (λ 2 = 0.68 µm) reflectances.In this case, the AART model gives an analytical expression for the optical snow grain diameter (Zege et al., 2008): Compared to the single-band algorithms for snow grain size retrievals (e.g.Nolin and Dozier, 1993;Stamnes et al., 2007;Tedesco and Kokhanovsky, 2007), which depend directly on the function R 0 and on its view-geometry-dependent errors, the band ratio algorithm to some extent suppresses these errors, reducing them to a second-order effect through the function A. We have done analysis with all of the previously described models (Sects.7.1-7.2) and found that the synthetic model with (σ = 0.3, D = 0.08) provides about the smallest  dispersion of retrieved grain size as well.Figure 9 shows the grain size retrieval with this model for experiment 2001b.
There is no regular dependence of retrieved grain size at relatively low view zenith angle θ < 35-40 • .At higher view angles (≥45-50 • ), the retrieved d systematically grows with θ 0 due to uncompensated BRF.Such dependence was observed from MODIS over Greenland (Lyapustin et al., 2009).From the reduced-scale image of Fig. 9, the average grain size at θ < 50 • can be assessed as ∼0.24±0.03mm.Further, the value d = 0.24 mm was used to study the spectral accuracy of the AART model.The results of this analysis are shown in Fig. 10.The three numbers above each image give the wavelength and parameters σ and D. In the low absorption region, λ ≤ 1.03 µm, the results are shown for the synthetic model with σ = 0.3 and D = 0.08.At medium absorption wavelengths, λ = 1.22 and 1.27 µm, a better fit was obtained with parameters σ = 0.2 and D = 0.10.At these wavelengths, the photon penetration depth is only a few millimeters.Thus, increasing absorption essentially increases the weight of the small-scale surface roughness (several mm vs several tens of cm) in averaging over the slope distribution function.
The maximal size of the high slope facets in the area of field measurements was only 5-10 cm (see also Fig. 1).With the average photon penetration depth of 30 cm in the visible spectrum, the assumption of independently scattering facets may not be entirely valid.Thus, the obtained value σ = 0.3 may be considered as a model result providing a good fit to the measurements.On the other hand, conditions for averaging are met in the near infrared, and value σ ∼ 0.2 may be better related to actual surface roughness.
At the high absorption wavelengths (λ = 1.65, 2.2 µm), the plane-parallel solution (σ = 0) with the synthetic model provides the best fit to the CAR BRF.In all wavelengths, the accuracy of the theoretical solution is close to the CAR BRF to within about ±0.05 over the full range of MODIS view angles θ ≤ 55-60 • , including the forward scattering and backscattering directions.

Conclusions
The spring 2008 ARCTAS experiment was one of major intensive field campaigns of the International Polar Year aimed at detailed characterization of atmospheric and surface properties of the Arctic region.In this work, we focused on processing and analysis of CAR-AATS-AERONET data from rather unique snow bidirectional reflectance P-3B airborne experiment.The goals of this study included obtaining snow bidirectional reflectance at high 1 • angular resolution from CAR measurements and using these data for an accuracy analysis of analytical RTLS, MRPV and AART BRF models over snow.Another major goal was developing a model of macroscopic surface roughness that would adjust the plane-parallel radiative transfer solution to experimental snow BRF.The main results of this work may be summarized as follows: 1.A method of atmospheric correction for CAR data was developed and thoroughly evaluated.The BRF and albedo are produced for 3 flight segments on 7 April and a flight segment on 15 April.The BRF and albedo are consistent on these two dates, and between different flight altitudes on the first date.Except for the red band, the derived CAR albedo is consistent with ground albedo measurements.The files containing the CAR measurements and derived spectral snow BRF can be found at http://car.gsfc.nasa.gov/data/index.php?mis id=8\&n=ARCTAS\&l=h.
2. The surface albedo derived from CAR data generally agreed well with ground measurements.
3. The obtained snow BRF is significantly anisotropic, even in the cross-plane.The derived BRF pattern from CAR measurements is similar to the measurements of Hudson et al. (2006) made over Antarctica from a 32 m tower.The use of Lambertian assumption for atmospheric correction may lead to large errors, especially in the shortwave channels.
4. Except for forward scattering (glint) region of angles ϕ < 40 • , the best fit MRPV and RTLS models provide a good fit to CAR BRF measurements to within ±0.05.
Overall, the RTLS model gives a better fit in the forward scattering angles, whereas the MRPV model suits snow reflectance better in the backscattering directions.
5. In agreement with the previous studies, the planeparallel radiative transfer solution was found to have large errors in the broad range of angles near the forward scattering and backscattering directions.Regardless of the shape of snow grains, the plane-parallel model significantly overestimates snow BRF in the broad glint region and underestimates it in the backscattering domain.Fitting CAR BRF data shows that the randomly oriented fractal model and the model of spheroids work much better than the Mie solution, and have complementary abilities in fitting the forward-and back-scattering regions.Based on this observation, we introduced a synthetic scattering function that is essentially a model of spheroids in the forward scattering angles, and is a fractal model in the backscattering range.The radiative transfer solution with synthetic phase function was found to fit measured BRF in the full range of angles better than any individual model.

6.
For the first time, we introduced averaging of the planeparallel radiative transfer solution over the slope distribution function that accounts for a natural snow surface roughness.Due to large snow grain sizes (compared to the wavelength), the scattering function of snow is rather flat in the backscattering domain and cannot provide the increase of reflectance at backscattering angles required to match observations.In these conditions, introducing rather steep slopes is perhaps the only way to increase snow reflectance at backscattering angles (see also Hudson and Warren, 2007).We found that averaging over slope distribution strongly reduces the difference between theoretical model and observations and allows us to model both the forward and backscattering regions with good accuracy at relatively high zenith angles.
7. Adding shadows even via a very simple model was found to further reduce the difference between the CAR data and the model.
8. A spectral-angular analysis showed that the AART model with the fitted surface roughness parameters σ and D provides an accuracy of about ±0.05 with the possible bias of ±0.03 in the spectral range 0.4-2.2µm at θ ≤ 55-60 • , including forward-and backscattering domains.
In future, the developed model of snow reflectance with surface roughness will be used to reduce BRF effect and improve snow grain size retrievals from MODIS with the band ratio technique (Lyapustin et al., 2009).
For the lack of better assumption, we will use an azimuthally independent Gaussian probability density function of slope distributions (Nakajima and Tanaka, 1983).
where µ n is cosine of zenith angle of the vector of normal.This model may not be valid in the world regions as Antarctica with predominant wind direction that creates sastrugi oriented perpendicular to the blowing wind (Hudson et al., 2006).Numerically, Eq. ( B1) is interpreted as follows.Let us define a reference right-handed coordinate system (x,y,z).Let the solar plane coincide with x-axes, ϕ 0 = 180 • .The vectors of incidence and reflection are defined by the directional cosines (or projections of a unit vector on each axis): for example, incidence vector is I = − 1 − µ 2 0 , 0, − µ 0 and view vector is V = − 1 − µ 2 cosϕ, 1 − µ 2 sinϕ, µ .
The minus sign at the x-projection accounts for the fact that the relative azimuth is defined with respect to the solar azimuth (ϕ = ϕ − ϕ 0 ).Next, for every orientation of a facet we define a new coordinate system (x n ,y n ,z n ) rotated in azimuth ϕ n about axis z and then rotated in angle θ n (µ n = cosθ n ) about axis y.The new coordinates of incidence and view vectors are related to the initial values by a rotation transformation, e.g.V 1 = T θ n T ϕ n V , where rotation matrices are To reduce computations and avoid numerical uncertainties caused by a definition of inverse cosine and inverse tangent, only new zenith angles (µ 01 , µ 1 ) are computed using rotation matrices while the new relative azimuth (ϕ 1 ) is calculated based on invariance of the scattering angle with respect to the rotation of the coordinate system.
Because zenith view and sun angles cannot be larger than 90 • in the rotated coordinate system (for plane-parallel model to hold), not every pair of zenith and azimuthal quadrature angles can be realized in Eq. (A1).In other words, for all geometries except nadir view and sun in zenith, only part of the upper hemisphere can be both illuminated and visible.This restriction is implemented through a joint constraint µ 01 ≤ 0, µ 1 ≥ 0. For this reason, the integrated solution needs a normalization coefficient, N = 2π 0 1 0 P (µ n ,ϕ n )dµ n dϕ n computed with the same constraint, which effectively shows the part of the hemisphere that is simultaneously illuminated and visible to the observer.The integration of Eq. ( B1) is performed numerically using Gaussian quadrature and a look-up table of function R precomputed with small step for the full range of azimuths and a range of zenith angles µ 0 , µ = 1 ÷ 0.12.The specific reflectance at quadrature points is obtained by 3-D interpolation in angles.

Shadow factor
The macroscopic surface roughness creates shadows that change bidirectional reflectance of snow.Even though shadowed areas are illuminated by the diffuse sky light, for simplicity we assume that shadow do not contribute into reflected signal.Here, we also ignore 3-D effects of horizontal diffusion of photons inside snow from illuminated into the shadowed area.This effect may be important in the visible part of spectrum, but it becomes negligible in the near infrared because of increased snow absorption.In this case, the measured signal will be reduced proportionally to the relative shadow area in the footprint: F Sh = (1 − S Sh /S).To compute shadow area, let us use a simple model of surface roughness as hemispheric protrusions of radius (height) r equidistantly spaced at a distance R. To simplify our model further, we assume that all of the shadow comes from illuminated half of the hemisphere (as if the sun were at horizon, θ 0 = 90 • ).In this case, the area of a shadow viewed at nadir is S Sh = π r 2 2 tgθ 0 .When the shadow is not subtended by the protrusion, the shadow area does not depend on the relative azimuth.This is approximately true for the forward scattering directions (|ϕ| ≤ 90 • ).In the backscattering directions, the shadow area will be reduced as S Sh ∼ π r 2 2 (tgθ 0 + tgθ cosϕ), where tgθ 0 + tgθ cosϕ ≥ 0. This formula is exact in the principal plane in the backscattering direction (ϕ = 180 • ) and in the cross-plane.Assuming that this expression is approximately valid for the other backscattering angles, we can write the final formula for the shadow factor: θ,ϕ), where H (θ 0 ,θ,ϕ) and D = r/R is density of protrusions.Although oversimplified, this shadow model predicts an overall correct angular dependence added by snow shadows.It has an advantage of dependence on a single parameter (density) that makes it rather robust for the remote sensing of snow given all sources of uncertainties and variability in modeling roughness of natural snow.Because of simplifications, such as neglecting diffuse irradiance of the shadowed areas, the best-fit parameter D matching experimentally measured BRF may differ from the real roughness density characterizing field conditions.The total snow reflectance accounting for slope distribution and shadows becomes: ρ(µ 0 ,µ,ϕ) = R new (µ 0 ,µ,ϕ)F Sh . (C2)

Figure 1 .Figure 2 .Fig. 1 .
Figure 1.Spectral albedo measurements on Elson Lagoon on April 15, 2008, with the P-3B airplane in the background.The leveled cosine collector is rotated to collect light from upward and downward hemispheres sequentially.The light is carried through a fiber-optic guide to the ASD spectral radiometer and its computer, which are mounted on the sled.

Figure 2 .Fig. 2 .
Figure 2. Comparison of the envelope of 15 April ground-measured albedos covering five site (lines) with CAR albedos obtained on 7 April (triangles) and 15 April (circles) at the Elson La goon site.

Fig. 3 .
Fig. 3. Snow BRF in CAR red band derived for the 7 April 2008 measurements at three different flight altitudes, and on 15 April at one elevation (d).
A.Lyapustin et al.:  Snow bidirectional reflectance from ARCTAS Spring-, and on 15 April at one elevation (d).

Fig. 5 .
Fig. 5. Retrieved snow BRF and accuracy of RTLS and MRPV models at 0.67 µm and 1.22 µm wavelength, for 7 April flight over Elson Lagoon.The accuracy, computed as (ρ − ρ Model ), is shown for the full range (on the left) and for a reduced range (±0.05, on the right) of values.

Fig. 6 .
Fig. 6.Error of atmospheric correction over snow using the Lambertian approximation (ρ − ρ Lamb ) for surface BRF, in visible and NIR wavelengths.

Fig
Fig. 7a-f.Difference between derived CAR BRF in the red band and modeled BRF (R new ) for different slope distribution functions.All images show differences with CAR red band BRF for 7 April flight over Elson Lagoon.

Fig. 7g -
Fig. 7g-m.Difference between derived CAR BRF in the red band and modeled BRF (R new ) for different slope distributions and sastrugi densities.The images (g-l) show results for flight 2001b (7 April), and the last image (m) shows result for flight 2005a (15 April).

Figure 7 Figure 8 .Fig. 8 .
Figure 7(g-m).Difference between derived CAR BRF in the red band and modeled BRF (R New ) for different slope distributions and sastrugi densities.The images (g-l) show results for flight 2001b (April 7), and the last image (m) shows result for flight 2005a (April 15).

Fig. 9 .
Fig. 9. Retrieved snow grain diameter from CAR BRF data for experiment 2001b.The right image shows the same data at reduced scale.

Fig. 10 .
Fig. 10.Difference between derived CAR BRF and modeled BRF (RNew) in seven CAR bands for flight 2001b.Three numbers above each image give the wavelength and parameters σ and D.

Table 1 .
Center wavelengths of CAR bands 3-8, 10, 13, in-band solar irradiance used in CAR calibration, solar irradiance based on the Spectral Irradiance Monitor (SIM) instrument on SORCE for 8 April 2008, and absorption optical thickness of atmospheric gases and ozone.

Table 2 .
Summary of CAR P-3B snow BRF experiments over Barrow.The last four columns give the average flight altitude, solar zenith angle and aerosol optical thickness of atmospheric column (AERONET) and of the layer above airplane (AATS) in the blue band.