Parameterizing anisotropic reflectance of snow surfaces from airborne digital camera observations in Antarctica

The surface reflection of solar radiation comprises an important boundary condition for solar radiative transfer simulations. In polar regions above snow surfaces, the surface reflection is particularly anisotropic due to low sun elevations and the highly anisotropic scattering phase function of the snow crystals. The characterization of this surface reflection anisotropy is essential for satellite remote sensing over both the Arctic and Antarctica. To quantify the angular snow reflection properties, the hemispherical-directional reflectance factor (HDRF) of snow surfaces was derived from airborne measurements in Antarctica 5 during austral summer in 2013/14. For this purpose, a digital 180° fish-eye camera (green channel, 490-585 nm wavelength band) was used. The HDRF was measured for different surface roughness conditions, optical-equivalent snow grain sizes, and solar zenith angles. The airborne observations covered an area of around 1000×1000 km in the vicinity of Kohnen station (75◦0′ S, 0◦4′ E) at the outer part of the East Antarctic Plateau. The observations include regions with higher (coastal areas) and lower (inner Antarctica) precipitation amounts and frequencies. The digital camera provided upward, angular10 dependent radiance measurements from the lower hemisphere. The comparison of the measured HDRF derived for smooth and rough snow surfaces (sastrugi) showed significant differences, which are superimposed on the diurnal cycle. By inverting a semi-empirical, kernel-driven bidirectional reflectance distribution function (BRDF) model, the measured HDRF of snow surfaces was parameterized as a function of solar zenith angle, surface roughness, and optical-equivalent snow grain size. This allows a direct comparison of the HDRF measurements with the BRDF derived from the Moderate Resolution Imaging 15 Spectroradiometer (MODIS) satellite product MCD43. For the analyzed cases, MODIS observations (545-565 nm wavelength band) generally underestimated the anisotropy of the surface reflection. The largest deviations were found for the volumetric model weight fvol (average underestimation by factor of 10). These deviations are likely linked to short-term changes in snow properties.

for accurately predicting future climate change in the polar regions. However, their remote and harsh environment makes them difficult to access and requires remote sensing techniques to observe snow surface reflection and its influencing parameters.
The high spatial and temporal variability of the snow surface albedo necessitates continuous observations with global coverage, which are provided by satellite instruments. Polar-orbiting satellites monitor the reflectance (i.e., reflected radiance in units of W m −2 sr −1 ) at the top of the atmosphere (TOA). However, they are restricted in terms of the number of available 5 observation angles and spectral bands as well as their temporal resolution. The processing of measurements from polar-orbiting satellites for monitoring the broadband surface albedo typically requires three steps (e.g., Qu et al., 2015): Atmospheric correction, modeling of the angular reflectance, and narrow-to-broadband conversion. During the first step, the TOA reflectance is converted into a surface reflectance by correcting for gaseous and aerosol scattering and absorption applying radiative transfer modeling (e.g., Vermote and Kotchenova, 2008). For the calculation of the narrowband surface albedo, an accurate knowledge 10 of the bidirectional reflectance distribution function (BRDF) of the surface is required due to the limited number of available observation angles. The BRDF describes the scattering of electromagnetic (EM) radiation from one incident direction into another direction of the hemisphere. By using a linear combination of the discrete spectral band measurements with specific weighting, the broadband surface albedo is calculated (e.g., Brest and Goward, 1987;Klein and Stroeve, 2002;Liang et al., 2003;Pohl et al., 2020). The largest uncertainty in this three-step process is introduced by the angular dependence of the 15 surface BRDF, especially when the reflection of the underlying surface is highly anisotropic. This is the case for the reflection of solar radiation at snow surfaces in polar regions due to low sun elevations and the anisotropic scattering phase function of snow crystals. As the snow surface of the Antarctic ice sheet is also used for the validation and cross-calibration of polarorbiting satellites (Jaross and Warner, 2008), a thorough understanding of the anisotropic reflection is needed. However, as an infinitesimal quantity, the BRDF of a surface cannot directly be measured under natural illumination conditions. From a 20 rigorous physical point of view, most satellite, airborne, and ground-based instruments measure the hemispherical-directional reflectance factor (HDRF, Schaepman-Strub et al., 2006) if the reflectance is constant over the full cone angle of the instrument's field of view (FOV). In contrast to the BRDF, the HDRF includes incident irradiance from the entire hemisphere.
An effective BRDF can only be derived from HDRF observations if the atmospheric influence is considered and the FOV is sufficiently small (e.g., Gatebe et al., 2003;Lyapustin et al., 2010). 25 For example, the retrieval of the aerosol optical depth (AOD) over highly reflecting surfaces is complicated, as it is hard to separate the surface and aerosol contributions to the measured reflectance at the TOA (Tomasi et al., 2015), especially for the comparably low values of AOD in polar regions (Mei et al., 2013). Due to its high surface albedo, most of the solar radiation incident on a snow surface is reflected, which depends strongly on the BRDF of the snow surface and, thus, on the illumination and viewing geometry. In general, most of the photons are scattered by the snow grains into the forward scattering direction. 30 Yang et al. (2014) compared Moderate Resolution Imaging Spectroradiometer (MODIS) measurements of AOD at 550 nm over Eastern China with ground-based data from the Aerosol Robotic Network (AERONET). Using a Lambertian model, 54 % of the regression points fell within the expected error envelope around the 1:1 line. This result is improved to 69 % if a more sophisticated, non-Lambertian model accounting for anisotropic reflection at the land surface is applied.
Yet, the BRDF varies with snow grain morphology, the solar zenith and azimuth angles, the liquid water content of the snowpack, the snow impurity type and concentration, the dimension and orientation of surface roughness structures, and the wavelength. To account for these effects, an accurate model of the snow BRDF is needed. A BRDF model can be either physical, empirical, or semi-empirical. Physical BRDF models (e.g., Cook and Torrance, 1982) accurately simulate the scattering of an EM wave at a surface by applying physical laws. The high accuracy comes at the cost of very complex computations. Empirical 5 BRDF models (e.g., Phong, 1975;Walthall et al., 1985) mimic the surface reflection by means of a simple, non-physical formulation. However, the drawback of the rather simple computations is their restricted accuracy. Semi-empirical models (e.g., Martonchik et al., 1998;Lucht et al., 2000) use simple, direct parameterizations of a more complex physical BRDF with a limited number of independent parameters.
The anisotropic reflection by snow surfaces was investigated by Li (1982) with simulations of the snow BRDF using Mie 10 theory and the doubling and adding method. Applying Mie theory and the discrete ordinate method, Han (1996) retrieved the surface albedo from satellite measurements in the Arctic.  developed a snow BRDF model including the effect of sastrugi by means of regularly spaced, identical, and rectangular protrusions.  and Leroux et al. (1999) employed the doubling and adding method, Mie theory, and ray tracing to develop a snow BRDF model including polarization. Comparing the simulated values with observations in the principal plane, they found that the BRDF in the near-15 infrared wavelength range is strongly affected by the snow grain shape, whereby simulations assuming hexagonal particle shapes yield an improved agreement with the observations compared to assuming spheres. Accordingly, Aoki et al. (2000) stressed the importance of the particle shape assumption for the scattering phase function used in their BRDF model (applying both Mie theory and the Henyey-Greenstein phase function) when comparing to observations between 0.52 µm and 2.21 µm wavelength. Kokhanovsky and Zege (2004) demonstrated the use of an asymptotic analytical equation to model the BRDF of 20 snow. Their approach represents the snow grains for example as fractal particles and, thus, accounts for their non-sphericity. Evaluating this asymptotic model with in situ measurements of the snow reflectance, Kokhanovsky et al. (2005) found a generally good agreement but a reduced model accuracy in the solar principal plane at large observation angles. In contrast to considering snow grains as independent scatterers of fractal shape, Malinka (2014) provided a framework to calculate the inherent optical properties (extinction coefficient, single scattering albedo, scattering phase function, and the polarization 25 properties) applying the concept of stereology that considers snow as a random mixture of ice and air. Malinka et al. (2016) combined the stereological approach with analytical, asymptotic formulas to calculate the bidirectional reflectance. Their model showed a high accuracy when compared to albedo observations of snow-covered sea ice.
The comparison of in situ measurements of the snow reflectance with simulations is essential in terms of model validation. Observations of the HDRF (or effective BRDF in case the FOV is small and the atmospheric influence is considered) 30 are conducted using a variety of different measurement concepts. For ground-based applications, manual or automated goniospectrometer systems are employed (e.g., Painter et al., 2003;Pegrum et al., 2006;Bourgeois et al., 2006a). Kuhn (1985) observed a peak in reflectance that is contained within ±60 • to both sides of the solar azimuth angle. This peak becomes more prominent with increasing solar zenith angle and snow grain size. Marks et al. (2015) measured the HDRF of snow surfaces between 400 and 1600 nm wavelength with the Gonio Radiometric Spectrometer System (GRASS, Pegrum et al.,35 2006) at Dome C, Antarctica. Their observations also showed enhanced forward scattering. In addition, they observed a larger anisotropy of the surface reflection at longer wavelengths. Employing a gonio-spectrometer, Bourgeois et al. (2006b) measured strong variations of the HDRF between 0.6 and 13 (in the wavelength range 350-1050 nm) depending on the solar zenith angle and the surface roughness at Summit, Greenland. Measurements with the Automated Spectro-Goniometer (ASG, Painter et al., 2003) showed a decrease of snow HDRF at all wavelengths between 400 to 2500 nm when the snow grain size increased 5 from 80 to 280 µm (Painter and Dozier, 2004b). Further comparisons with the results of a forward discrete ordinates radiative transfer model (Stamnes et al., 1988) revealed larger deviations between the simulations and observations for more complex crystals. The importance of the crystal habit for the anisotropy of reflection at snow surfaces was also emphasized by Dumont et al. (2010) and Stanton et al. (2016). The latter measured increasing anisotropy of the surface reflection during the growth of surface hoar in the laboratory. Several studies observed systematically less anisotropy for a typical snow HDRF than estimated 10 from simulations (Warren et al., 1998;Painter and Dozier, 2004a;Hudson et al., 2006;Hudson and Warren, 2007). In the solar principal plane, the models mainly overestimate the forward scattering and underestimate the backward scattering. Implementing non-spherical grains in the BRDF models (e.g., Mishchenko et al., 1999;Xie et al., 2006;Jin et al., 2008) improves the comparison with observations. The nonspherical model of Jin et al. (2008) agrees within ± 10 % for viewing zenith angles less than 60 • with observations in Antarctica performed by Hudson et al. (2006). However, the asymmetry between forward and 15 backward scattering still remains. This highlights the need to further incorporate macroscopic effects such as the roughness of the snow surfaces into the models Hudson and Warren, 2007).
Ground-based instruments observe the directional reflectance of a characteristic, homogenous surface, whereas airborne and satellite observations average over a larger measurement area. Thus, the latter are more suitable for studying the influence of macroscopic surface roughness on the surface HDRF. Nolin et al. (2002) employed multiangular measurements with the 20 multi-angle imaging spectroradiometer (MISR) for the characterization of surface roughness over Greenland and Antarctica.
Measurements with the clouds and the Earth's radiant energy system (CERES) showed monthly regional biases between -12 and 7.5 W m −2 in the cloudless TOA solar irradiance. These biases were attributed to the effect of sastrugi, which introduce a significant solar azimuth dependence: Kuchiki et al. (2011) observed a diurnal cycle in MODIS reflectances over the South Pole. In general, sastrugi decrease the forward scattering by casting shadows, increase the backward scattering due to a lower 25 effective incident angle caused by the sastrugi slope, and the snow HDRF loses its azimuthal symmetry (Warren et al., 1998). Zhuravleva and Kokhanovsky (2011) simulated a larger effect for a higher density of the sastrugi field.
Airborne measurements can be employed for the validation of BRDF models and the comparison with the large pixel size of satellite observations. Gatebe and King (2016) provided an extensive database of airborne spectral (effective) BRDFs for various surface types, e.g., ocean, vegetation, snow, desert, and clouds. The effective BRDFs were acquired by the Cloud 30 Absorption Radiometer (CAR, Gatebe et al., 2003) over a 30-year period between 1984 and 2014. The CAR is a scanning radiometer covering 14 spectral channels between 340 and 2324 nm. The effect of surface roughness on the effective BRDF was studied by Lyapustin et al. (2010) with CAR measurements during the Arctic Research of the Composition of the Troposphere from Aircraft and Satellites (ARCTAS) spring campaign in April 2008 (Jacob et al., 2010;Matsui et al., 2011). Their results showed an agreement between the kernel-based Ross-Thick Li-Sparse-Reciprocal (RTLSR) BRDF model (Lucht et al., 2000) 35 used in the operational MODIS BRDF/albedo product (Schaaf et al., 2002) and the CAR measurements to within ± 0.05.
The RTLSR model performed better in the forward scattering direction, whereas the kernel-based Modified Rahman-Pinty-Verstraete (MRPV) BRDF model (Martonchik et al., 1998) used for the processing of data from MISR performed better in the backscattering direction. Jiao et al. (2019) proposed an additional snow kernel within the kernel-driven BRDF model framework to better account for the anisotropic reflection of pure snow surfaces. Incorporating the additional snow kernel 5 yielded a correlation coefficient of 0.9 (compared to 0.65 for the original 3-kernel model) and only small biases between the model and different BRDF validation data from ground-based and satellite observations. Cox and Munk (1954) analyzed radiance calibrated analog photographs for the parameterization of the ocean HDRF. Nowadays, digital cameras are increasingly applied in vegetation and soil monitoring (e.g., Lebourgeois et al., 2008;Koukal and Atzberger, 2012). The instantaneous measurement of multiple viewing angles facilitates aerial HDRF measurements with dig- This study presents a method to derive the snow HDRF at visible wavelength (spectral band 490-585 nm) from airborne digital camera measurements in Antarctica during austral summer in 2013/14. Concurrent measurements of the optical-equivalent snow grain size retrieved from spectral surface albedo measurements, and surface roughness determined by means of a laser scanner, allow for investigating their effect on the snow HDRF in separate case studies. Subsequently, the measurements are 20 used to parameterize the snow HDRF applying the RTLSR model (Lucht et al., 2000). The presented methodological approach allows for a direct comparison of the airborne HDRF measurements with satellite observations from MODIS.
The definition of reflectance quantities used within this work, the modeling of the BRDF, and the inversion of a semiempirical, kernel-driven BRDF model are presented in Sect. 2. The field work and instrumentation are presented in Sect. 3 together with the detailed calibration involved in the measurement of the snow HDRF. The dependence of the snow HDRF 25 with respect to the solar zenith angle, the surface roughness, and the optical-equivalent snow grain size is studied in Sect. 4 based on three cases. Subsequently, the airborne measurements are compared with satellite-derived BRDF from MODIS. In the concluding Sect. 5, the findings of this work are summarized and perspectives for future studies are given.  The definition of the reflectance quantities applied within this work follows Schaepman-Strub et al. (2006). The intrinsic reflectance properties of a surface are given by its BRDF. It quantifies the reflection of the incident radiation at the surface and its scattering from one into another direction of the hemisphere. The spectral BRDF (f BRDF , unit of sr −1 ) provides for each zenith (θ i ) and azimuth angle (ϕ i ) of incident direct irradiance F i (θ i , ϕ i ; λ) the reflected radiance I r (θ i , ϕ i ; θ r , ϕ r ; λ) for all directions of reflection (defined by the reflection zenith and azimuth angles θ r and ϕ r ) by: The bidirectional reflectance factor (BRF, dimensionless) is obtained when the BRDF of a sample surface is divided by the BRDF of a diffuse (Lambertian) standard surface illuminated under the same conditions (identical beam geometry). As the 5 BRDF of an ideal Lambertian surface is (π sr) −1 , the BRF is given by In Eq. 2 and in the remainder of this section, the spectral dependence is omitted for reasons of simplicity. In atmospheric conditions, both the BRDF and BRF cannot be measured directly as the global irradiance F 0 reaching the surface is composed of a direct (F i ) and diffuse (F diff ) component. In this case, the measurable quantity is the HDRF. The definition of the HDRF 10 is analogous to the BRF, but includes irradiance from the entire hemisphere (denoted with 2π): (3) f dir denotes the fraction of direct incident radiation (i.e. f dir ∈ [0, 1]). The second step in Eq. 3 assumes that the incident diffuse radiation is isotropic (for details, see Schaepman-Strub et al., 2006). The additional integration over all reflection angles leads to the bihemispherical reflectance, generally called surface albedo α.

Modeling of the bidirectional reflectance
The shape of the BRDF is described by a weighted sum of trigonometric functions, generally referred to as kernels for volumetric scattering (K vol ), geometric scattering (K geo ), and isotropic scattering (K iso ).
The kernel-driven semi-empirical Ross-Li model (Lucht et al., 2000), which forms the basis of the MODIS 16-day BRDF/albedo product (Schaaf et al., 2002), is applied within this study. The BRDF is given as a linear combination of the 20 kernels with corresponding non-negative weighting functions f iso , f vol , and f geo : with the viewing zenith angle θ r , solar zenith angle θ 0 , relative viewing azimuth angle ∆ϕ, and wavelength λ. The kernels, depend on the scattering angle θ and the function O: cos θ = cos θ 0 · cos θ r + sin θ 0 · sin θ r · cos ∆ϕ, The functions C and D depend solely on the viewing and illumination geometry: Originally, the concept of the model was developed studying typical features in observations of the bidirectional reflectance of various surface types (Roujean et al., 1992). First, bare soil surfaces exhibit strong backscattering characteristics and show the effect of geometrical structures on the surface. Secondly, dense leaf canopies typically feature a minimum reflectance close to the nadir-direction that increases off-nadir for all azimuthal directions. In their approach, Roujean et al. (1992) represented 10 the BRDF model as a linear combination of these two observational characteristics; the former being accounted for by K geo , the latter by K vol . The volumetric kernel K vol stems from volume scattering radiative transfer models (Ross, 1981). Thereby, randomly located facets that absorb and scatter incident radiation are modeled under the single-scattering approximation. The formula given in Eq. 6 corresponds to the Ross-Thick kernel for a dense leaf canopy. The geometric kernel K geo is derived from surface scattering and geometric shadow casting theory (Li and Strahler, 1992) and expresses effects caused by intercrown

Inversion of semi-empirical, kernel-driven model
The main benefit of retrieving the weighting functions f iso , f vol , and f geo from the HDRF measurements is that they encompass all of the angular radiance information. The retrieval is done by inverting the RTLSR BRDF model. Thus, information about the complete two-dimensional shape of the HDRF is used.
The modeled f BRDF (θ i , θ r , ∆ϕ, λ) from Eq. 4 can be written in the form of a sum as: The spectral dependence is omitted here and the index l denotes a particular viewing and illumination geometry (θ i , θ r , ∆ϕ) l .
Considering an observation with N directional measurements ρ l (l = 1, . . . , N ), the error function E 2 is defined as the difference between the observed and the modeled reflectances such that, The degree of freedom d is (N − 3) and w l denotes weights which are assigned to the respective observations. In general, w l could take the values 1, ρ l , or ρ 2 l . The goal of the inversion is the determination of the model weighting functions f k such that E 2 is minimized. Strahler et al. (1996), amongst others, presented the analytical solutions for the three f k following this minimization and Lewis (1995) showed that E 2 has a global minimum in f k . 10 The inversion depends on the choice of the error function E 2 and, thus, on the choice of the weights w l . Minimizing the absolute error (using weights equal to unity) leads to smaller errors in angular domains with a large reflectance. In contrast, minimizing the relative error (e.g., w l = ρ l , ρ 2 l ) performs better in angular regions with lower reflectance. Subsequently, different choices for w l were tested for the retrievals performed in this work and w l = ρ 2 l was chosen as it produced the lowest retrieval errors across the entire angular domain.

Measurements and instrumentation
The airborne measurements were performed with the Polar 6 research aircraft (Wesche et al., 2016)

Surface roughness measurements using a laser scanner
The snow surface topography was measured using the airborne laser scanner RIEGL VQ580. In time-of-flight laser ranging, a near-infrared laser beam (1064 nm wavelength) is emitted downward and subsequently reflected upward by the snow surface before the echo is acquired by the sensor. From the time lag between emission and detection, the distance to the ground is calculated with a precision of about 25 mm depending on flight altitude. A fast-rotating polygonal mirror with a FOV of 60°a nd 10 to 150 scans per second allows for fully linear, unidirectional, and parallel scan lines.
After a correction for the aircraft attitude, a 1×1 km 2 digital elevation model (DEM) is generated from the geotagged laser point cloud with a resolution of 1 m. Subtracting the large scale topography (smoothed DEM), the residual field contains 5 the roughness information. The standard deviation of the residual field is interpreted as the surface roughness at the central coordinate of the DEM. Thus, roughness data is given for one data point per 1 km along the flight track of the aircraft. The uncertainty of the absolute height measurements (used for DEM generation) is less than 0.1 m. The relative analysis applied to the measurements yields an even higher accuracy.
The range in surface roughness is illustrated in Fig. 2b in the form of a histogram. The snow surface was generally smooth 10 with roughness structures mostly below 5-10 cm (note the logarithmic scale, maximum: 2.2 m).
3.2 Optical-equivalent snow grain size as retrieved from spectral surface albedo measurements Solar spectral radiation measurements were conducted using SMART (Wendisch et al., 2001;Ehrlich et al., 2008). The upward and downward spectral irradiance [F ↑ (λ), F ↓ (λ)] is measured within 0.3 to 2.2 µm wavelength applying an active horizontal adjustment system of the optical inlets to compensate for aircraft movement. From the irradiance measurements, the spectral 15 snow surface albedo was obtained (Wendisch et al., 2004). The spectral resolution is 2 to 3 nm between 0.3 and 1.0 µm and 15 nm up to 2.2 µm wavelength; the temporal resolution is in the order of 1 s. Combining the errors associated with the signal- to-noise ratio (1.3-3.0 %), the accuracy of the dark correction (0.1 %), the wavelength calibration (1.0 %), the accuracy of the cross-calibration (1.0-4.5 %), the non-ideal cosine response of the optical inlets (3.5 %), and the horizontal stabilization (1.0 %), the surface albedo measurements with SMART have an estimated uncertainty between 4.1 % and 8.1 % depending on wavelength (Carlsen et al., 2017).
The optical-equivalent snow grain size is defined as the radius R opt of a collection of spheres with the same volume-to- 5 surface ratio compared to the actual non-spherical snow grains. The retrieval of R opt from spectral surface albedo measurements is described in Carlsen et al. (2017). In principle, the snow grain size and pollution amount (SGSP) algorithm  was extended to spectral ratios of surface albedo at 1280 and 1100 nm wavelength. Being independent of systematic measurement uncertainties (e.g., cross-calibration of the optical inlets), this approach decreases the uncertainty of the retrieved R opt compared to the single-wavelength approach. The retrieved R opt from SMART and analogous ground-10 based measurements were compared to grain size observations utilizing reflectance measurements with MODIS (Carlsen et al., 2017). The grain size as retrieved from SMART measurements (mean value: 105 µm) is slightly higher than the MODIS retrievals (89 µm). The larger dataset of ground-based observations (72 µm) agrees well with the MODIS retrievals within the ranges given by the measurement uncertainties (linear correlation coefficient: 0.78, root-mean-square error: 24 µm).
The variability of R opt during the research flights is illustrated in Fig. 2c in the form of a histogram. The optical-equivalent snow grain size varied between 16 and 480 µm, but was mostly below 120 µm depending on precipitation and snow metamorphism processes.
3.3 Directional radiance measurements using a digital camera 3.3.1 Camera specifications 5 The digital camera Canon EOS-1D Mark III was used for the airborne HDRF measurements of the snow surface. It is a digital single-lens reflex camera with a complementary metal oxide semiconductor (CMOS) image sensor. The CMOS sensor covers 3908 × 2600 picture elements (pixels) on a sensor area of 28.1 × 18.7 mm (Advanced Photo System APS-H format). During the observations, the camera was configured with the 8 mm F3.5 EX DG Circular Fish-eye lens by Sigma. Color information is obtained by means of a color filter array (CFA) in front of the CMOS sensor. The CFA consists of red (R), green (G), and 10 blue (B) color filters, which determine the spectral response of the underlying pixel.
In principle, the measured signal s (in digital numbers, DN) should be saved essentially verbatim in the proprietary Raw format CR2 (Canon Raw version 2). This would ensure highest flexibility in the data analysis as no interpolation or white balance color correction is applied to the raw data. Full control over the post-processing steps is a prerequisite for radiometric measurements. During the airborne measurements, the camera was set to produce sRaw output format (small raw) with a 15 resolution of 1944 × 1296 pixels. At some point, internal chrominance subsampling is applied by the camera (Kerr, 2015).
The sRaw photos are decoded utilizing the open source tool dcraw (Coffin, 2017). The darkness level was set to 0 DN (see below) and the saturation level to 16383 DN as the camera captures images with color depth of 14 bit. In between, the raw data is linearly interpolated. The multipliers for all channels are set to one meaning no white balance color correction is applied.
The dynamic range of the output file is 16 bit (saturation at 65536 DN). 20 To characterize the linearity of the sensor, photos from the radiation emitted by an integrating sphere with varying output intensities were taken with the camera in the laboratory (Carlsen, 2018). The camera sensor showed a linear sensitivity over a large part of the dynamic range up to 53000 DN, the coefficient of correlation was exceptionally high with 0.99998. Photographs that measured higher signals (0.4 % of total number) were excluded from the analysis.
Compared to the dynamic range, the dark current (5 DN) and the read-out noise (7 DN) can be neglected. 25 Information on the spectral response of the camera is needed for the radiometric calibration. It determines the extent to which radiation of a certain wavelength passes the fish-eye lens as well as the CFA and is registered by the photodiodes. In this regard, the relative spectral response (RSR) function in units nm −1 of the three camera channels is defined as The RSR function of the camera was measured in the laboratory using an integrating sphere as a radiation source and varying the outgoing radiation between 300 and 750 nm in increments of 5 nm wavelength by means of a grating monochromator. The FOV of the camera is 180°due to the optics of the fish-eye lens. However, the camera is installed slightly above the lower aircraft body so that the camera is protected especially during take-off and landing. Therefore, parts of the aircraft frame 10 are constantly in the FOV of the camera, which is why the effective FOV is reduced to approximately 160°.

Radiance calibration and image post-processing
The post-processing of each raw image involves (a) the radiometric calibration, (b) the geometric calibration, (c) the aircraft attitude correction, and (d) the calculation of the HDRF. A detailed description of the different steps can be found in Carlsen (2018).

15
The measured signal was converted into the physical quantity of radiance (unit of W m −2 nm −1 sr −1 ) by means of a radiometric calibration in the laboratory. The pixel respond differently to an uniform illumination due to manufacturing tolerances (e.g., irregularities in the used silicon), contamination with dust particles, and optical effects at the edges of the lenses. These deviations lead to the photo response non-uniformity (PRNU) and need to be corrected.
The camera was positioned in front of an integrating sphere that served as an uniform radiation source. However, a distinct decrease in brightness is visible from the center to the edges of the sensor. This vignetting effect is typical for digital cameras (e.g., Lebourgeois et al., 2008). From the data sheet, the spectral radiance I sphere (λ) that the integrating sphere emits at the 5 specific optometer current is known. The calibration factor k c is defined at each pixel location on the sensor (row x, column y) and carries the unit of W m −2 nm −1 sr −1 (DN/s) −1 . During the calibration, the exposure time t exp was set to 1/1000 s.
Not only does k c correct for the PRNU, it simultaneously performs the absolute calibration transforming the measured digital 10 signal into the physical quantity of radiance with units W m −2 nm −1 sr −1 .
The geometric calibration of the camera-lens system relates each sensor pixel to a viewing zenith and azimuth angle (θ v , ϕ v ). Often, the process of geometric calibration of a camera involves calibration equipment or the use of planar targets such as checkerboard patterns (e.g., Tsai, 1987;Urquhart et al., 2016). Within this work, a stellar calibration method is applied (e.g., Schmid, 1974;Klaus et al., 2004;Mori et al., 2013;Urquhart et al., 2016) utilizing the high precision to which the positions 15 of stellar objects are known. In this regard, 684 different star positions identified in subsequent pictures of the night sky were utilized to calculate θ v and ϕ v for each sensor pixel.

Calculation and uncertainties of measured surface HDRF
The camera took pictures with a temporal resolution of about 8 s. Exemplarily, Fig. 4 demonstrates the derivation of the snow HDRF from a raw image taken by the camera on 2 January 2014 (8:16 UTC, see Fig. 4a). First, the pixels receiving 20 radiation from the direction of the aircraft frame are excluded. For each pixel location (x, y), the radiance I(x, y) (in units W m −2 nm −1 sr −1 ) is calculated from the measured signal s using the absolute calibration factor k c and the exposure time t exp : The camera viewing angles are calculated from the geometric calibration for each pixel (x, y). As the camera is fixed to the 25 aircraft frame, a correction for the aircraft attitude was implemented to obtain the reflection angles θ r and ϕ r . Thus, beside the geometric calibration, the observed reflection angles are determined by the attitude angles of the aircraft. Utilizing the data from the internal navigation system (INS) and the global positioning system (GPS) on Polar 6, the viewing angles are corrected depending on the roll, pitch and yaw angles of the aircraft at the time of measurement. In this regard, Euler rotations are applied as described in Ehrlich et al. (2012). Finally, the observed HDRF is calculated by: 30 R HDRF (θ 0 , ϕ 0 ; θ r , ϕ r ) = π sr · I(θ r , ϕ r ) F ↓ (θ 0 , ϕ 0 ) .
ϕ 0 ). Based on the resulting reflection angles, a polar plot of the measured R HDRF is created (camera channel G, see Fig. 4b).
To achieve comparability, each image is rotated into the azimuthal direction of the Sun.
The downward irradiance measurements from SMART could not be used for the calculation of the HDRF due to calibration issues. Note that this pertains to the radiometric calibration only, which converts the digital numbers registered by the 5 spectrometer into units of irradiance. For the albedo measurements with SMART, a relative calibration of the upper and lower sensors is sufficient and an absolute radiometric calibration is not required. Thus, the albedo measurements and the retrieval of the optical-equivalent snow grain size are unaffected by this calibration issues. For the calculation of the HDRF, the global irradiance was simulated along the flight track with the library for radiative transfer libRadtran (Mayer and Kylling, 2005) using the discrete ordinate radiative transfer solver DISORT by Stamnes et al. (1988) instead. Vertical profile data of air temperature the global irradiance F ↓ in Eq. (18). The uncertainties in the HDRF measurements stem from sensor characteristics (estimated with 0.5 % due to signal-to-noise ratio, dark current, linearity, read-out noise, chrominance subsampling), the radiometric calibration (4 %), the geometric calibration and the correction for the aircraft attitude (combined estimate of 1.0 %).

Averaging
From trigonometric considerations, the footprint of the camera (twice the radius r of the disc on the ground pictured by  (e) Relative difference (ρ l − R HDRF,l ) /R HDRF,l . The accumulated RMSE is 0.04.

Approximation of surface BRF with HDRF measurements
For each averaged HDRF, the weighting functions f iso , f vol , and f geo are retrieved by inverting the RTLSR BRDF model and setting the HDRF equal to the BRDF in Eq. 12. This is necessary as the BRDF is not measurable under atmospheric conditions. However, the atmospheric influence is assumed to be small due to the high surface elevation and the dry, aerosolfree atmosphere on the Antarctic plateau: the mean AOD at Kohnen station was 0.015 between 2001 and 2006 (at 500 nm, simulated the difference between the HDRF and BRF (R BRF = π · f BRDF , compare Eq. 2) for snow surfaces for different fractions f dir at 550 nm wavelength. They found that with an increasing diffuse fraction of the incident irradiance, the shape of the HDRF is smoothed in comparison to the BRF (f dir = 1, see Eq. 3). For f dir = 0.8 at θ 0 = 30 • , the shape of the HDRF is still close to that of the BRF. Hence, the measured HDRF shapes can serve as approximations for the intrinsic BRF of the 15 underlying snow surface for the measurements analyzed within this study (f dir > 0.8). (see Fig. 4e). The largest deviations occur in the forward scattering direction (up to 3 %) because the location of maximum is not mimicked perfectly by the modeled HDRF. However, up to viewing zenith angles of around 75 • , the relative difference is below 1 %.

Quality of the inversion
The root-mean-square error (RMSE), 5 serves as a criterion for the quality of the inversion. Only RMSE below 0.1 are considered full inversions within the MODIS retrieval (e.g., Stroeve et al., 2005). For the example shown in Fig. 4, it is 4.0 %. σ RMSE strongly depends on the anisotropy of the observed HDRF. Figure 5 shows the RMSE for all airborne retrievals as a histogram, with 96 % of the retrievals showing an RMSE below 0.1. Hence, the upper bound for the uncertainty of the inversion is in the range of 10 %.
The influence of the angular sampling on the retrieval of the model parameters can be investigated using the weights of  f k is negligible. The mean standard deviations of 0.031 (f iso ), 0.020 (f vol ), and 0.003 (f geo ) are used as error bars for the f k within this study.

Comparison with BRDF derived from MODIS satellite observations
The use of the RTLSR model in the retrieval of the model parameters f k from airborne observations provides the framework to compare with satellite observations. The MODIS BRDF/Albedo Model Parameters product (MCD43A1) provides the f k in 500 m resolution utilizing multi-date, atmospherically corrected, and cloud-cleared input data over a period of 16 days (Schaaf et al., 2002). Version 6 of the MCD43A1 product uses both Terra and Aqua data and is temporally weighted to the 9th day of a 5 16-day retrieval window, thereby putting higher emphasis on the actual day of interest than previous versions. The combination of MODIS on Terra and Aqua provides a higher number of high quality reflectance measurements, resulting in better temporal and angular sampling. As the retrieval uncertainty increases with larger solar zenith angle, a thorough quality assessment of the satellite data is needed especially in polar regions. Stroeve et al. (2005) compared the MODIS albedo product with in situ observations on the Greenland ice sheet and recommended to only use full inversion data (RMSE < 0.1, WoD < 2.5). Hence, 10 this study restricts the satellite data to full inversions based on the quality flag information provided in the MCD43A2 product (quality flag values 0 or 1). Lower quality magnitude inversions would facilitate archetypal BRDF parameters to improve the retrieval in case of less than 7 high quality reflectance observations or poor angular sampling. The BRDF data from MODIS spectral band 4 (0.545-0.565 µm) is used as it coincides best with the green camera channel (see Fig. 3).
The MODIS data was resampled on the flight track using nearest neighbour with respect to the great circle distance. For

Results and discussion
To study the influence of the solar zenith angle, the surface roughness, and the optical-equivalent snow grain size on the snow HDRF, the effects of the individual parameters need to be separated. For this purpose, several case studies were selected: (a) measurements at different θ 0 (marked with red circles in Fig. 2a

Influence of solar zenith angle
The research flights with Polar 6 on 28 December 2013 were exploited to study the influence of θ 0 on the snow HDRF. Five consecutive flights were conducted in an area northeast of Kohnen (see red flight track in Fig. 1a).
For each flight, the snow HDRF was averaged over a 30 s-segment during the easternmost flight leg; the results are shown in the top panel of Fig. 6. The solar zenith angle varied between 51.8°and 71.6°. As Polar 6 took identical routes throughout the 5 flight, the respective flight leg always covered the same area. The central points of the five consecutive segments are separated by less than 1 km. Therefore, the surface roughness (retrieved from the laser scanner) of approximately 6 to 7 cm remained constant. Within the retrieval uncertainties, the same holds true for the optical-equivalent snow grain size. R opt varied between 70 µm and 85 µm. This minimizes a possible influence of R opt or l rough on the snow HDRF that would be superimposed on the effect of θ 0 . Cloudless conditions prevailed during the flights. 10 With increasing θ 0 , the maximum of the HDRF in the forward scattering direction becomes more pronounced (see Fig. 6).
Correspondingly, the anisotropy gets larger. Figure 6 shows the dependence of the f k from θ 0 for the green camera channel (490-585 nm wavelength). f iso shows no clear dependence on θ 0 and varies between 1.06 and 1.10. f geo weakly increases with θ 0 from 0.03 at θ 0 = 51.8 • to 0.05 at 71.6°. During the first three research flights, θ 0 did vary only slightly. Therefore, the strongest trends are visible between flights 3 and 5.
In particular, f vol increases from 0.25 to 0.36.

15
The probability for photons entering the snowpack to leave it after just a few scattering events increases for lower Sun elevation. With increasing θ 0 , the reflection properties of the snow layer converge to the single-scattering properties of ice crystals.
In addition, longer shadows are cast by surface roughness structures at larger θ 0 . This expected increase in the anisotropy of the snow HDRF for increasing θ 0 is obvious in the changing model parameters f k and has been observed in earlier studies (e.g., Dirmhirn and Eaton, 1975;Kuhn, 1985;Warren et al., 1998;Hudson et al., 2006). For the smooth surfaces, the maximum in forward scattering direction is more pronounced. With increasing surface rough-30 ness, the anisotropy gets lower. Roughness structures enhance the backscatter by changing the effective angle of incidence. In addition, they cast shadows that reduce the forward scatter. The evolution of the retrieved model parameters f k as illustrated in Fig. 7 demonstrates the decreasing anisotropy: f iso and f geo show a slight decreasing trend whereas, most prominently, f vol decreases from 0.33 to 0.24 with increasing surface roughness. This is in accordance with Warren et al. (1998) who observed a reduction of the forward reflection peak due to sastrugi during tower measurements at South Pole Station. The effect of surface roughness on the BRDF has been previously investigated in observational (e.g., Grenfell et al., 1994;Hudson and Warren, 2007) and modeling studies (e.g., Zhuravleva and Kokhanovsky, 2011) and is facilitated by the multi-angular reflectance data from MISR to map surface roughness over ice sheets and sea ice from space (Nolin et al., 2002).

Influence of optical-equivalent snow grain size
The separation of the two parameters θ 0 and l rough showed that they influence the snow HDRF in opposite ways. Ice absorption becomes dominant in the near-infrared part of the solar EM spectrum. Hence, even though a wide variety of optical-equivalent snow grain sizes was covered during the flights (16 to 480 µm), no large effect on the snow HDRF was expected as the camera is only sensitive in the visible wavelength range (Warren et al., 1998, compare Fig. 3). However, 32 observations between 5 14:45 and 15:15 UTC on 2 January 2014 (see brown flight track in Fig. 1a) and depicted in Fig. 8 show a decrease of f vol from 0.34 at 184 µm to 0.17 at 325 µm snow grain size. This is consistent with the decrease in anisotropy as visible in the averaged HDRF measurements (top panel in Fig. 8).
However, this stands in contrast to findings by Steffen (1987), who measured the bidirectional reflectance on a glacier in Northwestern China for powder snow (grain size 0.15 mm), new snow (0.5-1.0 mm), and old snow (1-3 mm) and found higher 10 anisotropy for larger snow grains. It should be noted that the range of snow grain sizes covered within this study is significantly smaller and the footprint of the measurements is larger compared to the local study by Steffen (1987). Even though other parameters were kept constant during this case study (θ 0 : 51.3-53.4 • , l rough : 4-7 cm, footprint: 1.3-1.7 km), effects of other influencing parameters cannot be ruled out completely. For example, Steffen (1987) reported changing ice crystal shapes between the measurements and suggested a possible influence of dust particles within the snow. This could partly explain the 15 opposite conclusions drawn within this study.
The linear correlation coefficients (R 2 ) are 0.98 (for θ 0 ), 0.75 (l rough ), and 0.68 (R opt ). Note that these parameterizations are only valid for the wavelength range of about 490-585 nm (green camera channel) and the ranges of all three parameters as given within the descriptions of the case studies.   Thirdly, the airborne measurements still included varying footprint sizes (mostly between 1 and 3 km) compared to the 500 m resolution of the MODIS product. This was the result of trading off comparability with MODIS and data coverage.
Also, no correlation between the f k retrieved from the camera measurements and the footprint size was found (not shown here).
Lastly, comparing the airborne camera data with the MODIS product poses challenges that are inherent to the measurement 25 setup. The observation times of the MODIS product do not coincide with the times of the research flights. MODIS spectral band 4 was chosen as its RSR function agrees best with the green camera channel. As the HDRF is wavelength dependent, differences naturally occur due to different spectral coverage but should be small.
Thus, despite the thorough quality assessment analogously applied to both airborne and satellite data, a quantitative comparison between the two retrievals seems challenging. Nevertheless, the comparison highlights the difficulties to accurately

Conclusions
The HDRF of snow surfaces in the wavelength band 490-585 nm was derived from airborne measurements in Antarctica during austral summer in 2013/14. An airborne digital 180°fish-eye camera installed on a research aircraft was used for this purpose. The surface roughness was retrieved using an airborne laser scanner, while the optical-equivalent snow grain size was calculated from spectral surface albedo measurements. While based at Kohnen station (75 • 0 S, 0 • 4 E) at the outer part of the 5 East Antarctic Plateau, the airborne measurements with the Polar 6 research aircraft covered an area of around 1000×1000 km 2 in Dronning Maud Land. Thus, the snow HDRF was obtained for a variety of conditions with different solar zenith angle (θ 0 ), surface roughness (l rough ), and optical-equivalent snow grain size (R opt ).
The digital camera provides airborne radiance measurements with high angular resolution. The technical characterization of the digital camera in the laboratory revealed excellent linearity as well as negligible dark current and read-out noise. From 10 the measured angular radiances, the HDRF was calculated applying simulations of the global irradiance with the library for radiative transfer libRadtran using the discrete ordinate radiative transfer solver DISORT. The relative uncertainty of the HDRF measurements is estimated with 4.5 %. The footprint of the snow HDRF measurements analyzed within this study varies between 1 km to 4.5 km.
Three case studies were investigated to separate the effects of θ 0 , l rough , and R opt on the snow HDRF. With increasing solar zenith angle, the HDRF maximum in the forward scattering direction becomes more pronounced. Conversely, the backscatter is enhanced over rougher surfaces by changing the effective angle of incidence and by casting shadows. This was quantified by inverting the semi-empirical kernel-driven RTLSR model and parameterizing the snow HDRF with respect to θ 0 , l rough , and R opt (Eqs. [20][21][22]. The uncertainty of the inversion is estimated with 10 %. The increased anisotropy (with increasing θ 0 ) is mainly shown by the increase in the model parameter f vol . Vice versa, f vol decreases with increasing l rough . The snow 20 grain size reveals a similar effect as surface roughness structures in terms of a decrease of the anisotropy with increasing R opt , which differs from earlier findings (e.g., Steffen, 1987). Possible reasons are discussed (footprint size, ice crystal shape, contamination with dust particles). However, to yield stronger dependence on the optical-equivalent snow grain size, a camera sensitive to radiation in the near-infrared part of the EM spectrum needs to be employed.
Continuous satellite observations of the snow surface albedo with global coverage critically depend on the angular modeling 5 of the surface BRDF, which may introduce large errors especially in polar regions due to low sun elevations and the anisotropic scattering phase function of snow crystals. The RTLSR model was chosen for the inversion as it forms the basis of the MODIS 16-day BRDF/albedo product. Thus, the retrieved model parameters f k from the HDRF measurements allow for direct comparison with BRDF products from satellite remote sensing. Despite the similar quality assessment applied to airborne and satellite data, large deviations were found especially for the volumetric model weight f vol . The airborne values for f vol are larger than 10 the corresponding MODIS retrievals (by at least a factor of 10). Short-term changes in snow properties (precipitation, drifting snow), which are not captured by the 16-day MODIS retrieval, and different solar zenith angles at the time of measurements are discussed as likely reasons for the observed differences. Although the effects are presumably small, influences of a varying footprint size for the airborne observations and the different RSR functions of the green camera channel and MODIS spectral band 4 can further contribute to this discrepancies.

15
Generally, MODIS observations underestimated the anisotropy of the reflection by the snow surfaces. This has possible implications for satellite retrievals of surface albedo and atmospheric parameters over snow surfaces (e.g., AOD, cloud properties), which strongly depend on a correct angular modeling of the surface reflectance (e.g., Qu et al., 2015).
The analysis of the snow HDRF measurements can be readily applied to other airborne campaigns. For example, the same digital camera with a similar setup of the Polar 6 research aircraft was used (a) during the Arctic CLoud Observations Using Lucht, W., Schaaf, C. B., and Strahler, A. H.: An algorithm for the retrieval of albedo from space using semiempirical BRDF models, IEEE