On the retrieval of snow grain morphology, the accuracy of simulated reflectance over snow using airborne measurements in the Arctic

Accurate knowledge of the reflectance from snow/ice covered surface is of fundamental importance for the investigation and retrieval of snow parameters and atmospheric constituents. This is a prerequisite for identifying and quantifying changes in the environment and climate in Polar Regions. However, the current differences between simulated and measured reflectance in a coupled snow-atmosphere system, leads to systematic errors in the determination of the amount of trace gases, aerosol and cloud parameters from space based and airborne passive remote sensing observations. In this paper, we describe studies of the retrieval of snow grain morphologies, also called habits, and their use to determine reflectance and test the accuracy of our radiative transfer model simulations of reflectance by comparison with measurements. Firstly we report on a sensitivity study. This addresses the requirement for adequate a priori knowledge about snow models and ancillary information about the atmosphere; the objective being to minimize differences between measurements and simulation. For this aim we use the well-validated phenomenological radiative transfer model SCIATRAN. Secondly and more importantly, we present a novel two-stage snow grain morphology (i.e. size and shape of ice crystals in the snow) retrieval algorithm. We then describe the use of this new retrieval to estimate the most representative snow model, using different types of snow morphologies, for the airborne observation conditions, performed by NASA’s Cloud Absorption Radiometer (CAR). The results show that the retrieved ice crystal shapes are consistent with the expected snow morphology (estimated from temperature information) in the measurement area over Barrow/Utqiaġvik, Alaska in 2008. Thirdly, we present a comprehensive comparison of the simulated reflectance (using retrieved snow grain size and shape as well as independent atmospheric data) with that from airborne CAR measurements in the visible and NIR wavelength range. The results of this comparison are used to assess the quality and accuracy of the radiative transfer model in the simulation of the reflectance in a coupled snow-atmosphere system. Assuming that that snow layer consists of ice crystals with aggregates of 8 column ice habit having an effective radius ~ 98.83 μm, we find that for a surface covered by old snow, the Pearson correlation coefficient, R, between measurements and

retrieval algorithm are presented. In sect. 7, the results of the reflectance simulations are compared to CAR measurements.
Finally, conclusions are drawn in sect. 8. Appendix contains detailed description of the snow grain size and shape retrieval algorithm used in the study.

Theoretical background
To describe the directional signature of reflectance over different surface types, the Bidirectional Reflectance Distribution Function (BRDF) as defined by Nicodemus (1965), is the commonly used reflectance quantity. The term BRDF describes the reflection of incident solar radiation from one direction to another direction (Nicodemus 1965). The mathematical form of BRDF is expressed as (Nicodemus et al., 1977;Schaepman-Strub et al., 2006): where Lr is the reflected radiance, θ and φ are the zenith and azimuth angles, respectively. The subscript i corresponds to the incident and r to the reflected beams. E is the incident surface flux (irradiance) and λ is the wavelength. However, the BRDF is not a directly measurable quantity because of its being formulated as a ratio of infinitesimal quantities (Nicodemus et al., 1977;Schaepman-Strub et al., 2006). Nicodemus et al. (1977) provided an extensive description of reflectance terminologies and measurable quantities e.g. the Bidirectional Reflectance Factor (BRF), the Hemispherical Directional Reflectance Factor (HDRF), the Directional Hemispherical reflectance (DHR), etc. According to Nicodemus et al. (1977) and Schaepman-Strub et al. (2006) each of the terms is defined for the specific illumination and reflectance geometries for which, the reflectance properties are measured e.g. satellite, airborne or laboratory measurement conditions. Following the method of Gatebe and King (2016), the effective BRDF at a horizontal (flat) reference plane is defined as: where BRDF λ e is as an average of the BRDF over an appropriate area, angle and solid angle for specific observation geometry; ∆ ω i is a finite solid angle element. The validity of this approximation relies on the experimental evidence that the BRDF is not significantly influenced by the following effects: the finite intervals of area, angle, solid angle and the distribution function; sub-surface scattering; radiation parameters such as wavelength and polarization, fluorescence etc. (i.e. significant variations do not occur within small intervals, see Nicodemus et al., 1977;Gatebe and King, 2016). As a result, the BRDF λ e is determined by: BRDF λ e = L r e (θ i , θ r , ∆ φ) where L r e is the measured radiance, F 0 , λ is the solar irradiance incident at the top of atmosphere (TOA). Often, it is helpful to have a description of the difference between the measured surface reflectance and a Lambertian reflector; in such a case the equivalent Bidirectional Reflectance Factor BRF λ e , which is BRDF λ e multiplied by π is more representative.
To isolate the reflectance properties of surface and derive BRF λ e or BRDF λ e just above the surface, we need to apply atmospheric correction methods on measured radiance at TOA or flight altitude e.g. by using knowledge of the atmospheric scattering or absorption using RTMs. This removes from the measured radiance, the four atmospheric contributions from the atmosphere at TOA or flight altitude (Schaepman-Strub et al., 2006): i) the atmospheric path radiance, ii) the scattering by the atmosphere before the solar radiation has reached the surface, iii) the scattering by the atmosphere after being reflected by the surface, iv) the scattering by the atmosphere before and after reaching the surface. However, most of the atmospheric contributions in measurements close to the surface are negligible (except diffuse component number ii) and measured quantities represent the "at surface" radiance (Schaepman-Strub et al., 2006). Sensitivity studies have demonstrated that atmospheric contributions to the CAR channel observations range from 3 to 12% depending on wavelength (Soulen et al., 2000). Consequently, previous studies presented either the BRFs in a surface-atmosphere system at flight altitude without atmospheric correction (Soulen et al., 2000) or the BRFs right above the surface after atmospheric correction (Gatebe et al., 2005;Gatebe and King 2016).
In this study, to avoid uncertainties arising from different assumptions being part of the atmospheric correction methods, no such correction is applied to measured radiances L r , h at flight altitude h. Instead, the reflectance at flight altitude in the snow-atmosphere system is calculated by the equation: where L r , h is the measured radiance at flight altitude. All reflectance/ BRF λ e values at flight altitude, presented in this study are calculated using Eq. 4 and referred to as "reflectance". In the simulation of the reflectance in a coupled snow-atmosphere system, to account for atmospheric effects contribution properly, independent data about atmospheric parameters (Aerosol Optical Thickness (AOT) and gases absorption) at the time and close to the location of measurements are needed and taken from ground-based and space-borne measurements and applied to the simulation. More details of the atmospheric data are discussed in sect. 3 and 4. To estimate BRF λ e just above the surface, further atmospheric correction is needed. We assume at infrared wavelengths where atmospheric scattering is negligible, the reflectance at flight altitude is a good estimation of BRF λ e just above the surface.

Measurements
CAR is an airborne instrument, developed at NASA's Goddard Space Flight Center. It has been used during several field campaigns around the world since 1984 up to present to measure the single scattering albedo of clouds and the bidirectional reflectance of various surface types etc. For this study, we used CAR data from the ARCTAS campaign conducted at Elson Lagoon, near Barrow/Utqiaġvik, Alaska, in April 2008 (Lyapustin et al., 2010;Gatebe and King, 2016). The goal of ARCTAS was to study physical and chemical processes in the Arctic atmosphere and surface. Date, location, measurement geometry and available atmospheric parameters during the measurements used in this study are presented in Table 1. The unique design of CAR provides simultaneously both up-welling and down-welling radiances at 14 spectral bands located in the atmospheric window regions of UV, visible and near-infrared from 0.34 μm to 2.3 μm comprising important wavelengths relevant for remote sensing applications such as aerosol retrievals. Through a rotating scan mirror, the instrument provides viewing geometries suitable for measurements needed for BRF calculation. CAR collects data bya mirror rotating 360° in a plane perpendicular to the direction of flight through a 190° aperture that allows acquiring data from local zenith to nadir or horizon to horizon with an angular resolution of 1°. The high angular/spatial resolution of 1° in both viewing and azimuth angles allowed the anisotropy of the reflectance in the snow-atmosphere system to be estimated with high accuracy. The spatial resolution of CAR depends on the flight altitude e.g. 10 m 2 and 18 m 2 at nadir for 600 m and 1000 m flight altitude, respectively, which increases with the viewing zenith angle (VZA) e.  Fig. 1 and Fig. 2. In spite of the influence of the atmospheric scattering and absorption, the general features of the snow BRF are clearly observable. The latter comprise: i) the decrease of snow reflectance with increasing wavelength due to the increasing absorption by snow at longer wavelengths; ii) the increase of the snow BRF as a function of VZA and the strong forward scattering peak in the principal plane at large VZA; iii) the smaller angular variation of the BRF at cross plane compared to the principal plane. The reflectance values increase with altitude. The snow surface inhomogeneity decreases with increasing altitude due to the change of spatial resolution with altitude (Gatebe and King, 2016;Lyapustin et al., 2010). Accordingly, at poorer spatial resolution, spatial homogeneity are more efficiently averaged.
To account for aerosols, we use the Aerosol optical thickness, AOT data acquired by the nearby Aerosol Robotic Network (AERONET) sun-photometer at Barrow/Utqiaġvik during the CAR measurement time. AERONET is a globally distributed network and provides long-term and continuous ground-based measurements of the total column aerosol optical thickness derived from the attenuation of sun light and provided often at high temporal resolution of 15 minutes. AERONET AOT data are provided at 0.5 μm and 0.6 μm wavelengths. We use the Ångström exponent to calculate AOT values at the reference wavelength (0.55 μm) required for the SCIATRAN simulation. Table 1 shows the calculated AOT at 0.55 μm based on AERONET data for Barrow/Utqiaġvik at the closest time to the CAR airborne measurements.
To account for ozone absorption, we use knowledge of the ozone total column amount retrieved from the space borne measurements by using the University of Bremen weighting function DOAS (WFDOAS) algorithm version 4 (Weber et al., 2018). The WFDOAS data are selected using the criteria of having smallest temporal and spatial differences with CAR data.
For nitrogen dioxide, we use vertical column information from the SCIATRAN database obtained from a 2D chemical transport model developed at University of Bremen (Sinnhuber et al., 2009).
The derived AOT and trace vertical column have been used in the simulation of radiative transfer processes in the snowatmosphere system.

Simulations
SCIATRAN is a software package for radiative transfer modeling, developed at the Institute of Environmental Physics, University of Bremen (Rozanov et al., 2002;2014) and freely available at http://www.iup.uni-bremen.de/sciatran/. The SCIATRAN package has been used in a variety of remote sensing studies to simulate radiative transfer processes in the spectral range from the ultraviolet to the thermal infrared (0.18 μm -40 μm), assuming either a plane parallel or a spherical atmosphere (Rozanov et al., 2014).
To calculate reflectance values, SCIATRAN assumes that the snow is a layer with an optical thickness of 1000 and a geometrical thickness of 1 m composed of ice crystals of different morphologies and placed above a black surface. This assumption was successfully validated by Rozanov et al. (2014). The snow layer is assumed to be vertically and horizontally homogeneous and composed of a monodisperse population of ice crystals. The impact of impurities in the snow e.g. dust, black carbon etc. is neglected in this study. To simulate the radiative transfer through a snow layer, the single scattering properties of ice crystals including extinction and scattering efficiencies, single scattering albedo and phase functions need to be defined in SCIATRAN. All these parameters are dependent on the wavelength, size and shape of the particle (Leroux et al., 1999). Recently, a new data library of basic single scattering properties of ice crystal habits developed by Yang et al. (2013) has been incorporated in the SCIATRAN model (Pohl et al., Personal communication). This database comprises a full set of single scattering properties at wavelengths from the UV to the far IR for the following eight ice crystal morphologies: droxtal, column and hollow column, aggregate of eight columns, plates, small aggregate of five plates, large aggregate of ten plates, and hollow bullet rosettes. More detailed information about the ice crystal shapes and sizes can be found in Yang et al. (2013). In addition to the above-mentioned eight ice crystals, optical parameters for triadic Koch fractal (referred as fractal in this paper) particles are used as well (Macke et al., 1996;Rozanov et al., 2014). The fractal particle model uses regular tetrahedrons as its basic elements. In this study, the second generation fractals as described in Macke et al. (1996) and Rozanov et al. (2014) are utilized.
In SCIATRAN, the snow grains are specified by their single-scattering properties of sparsely distributed particles.
Namely, the snow grains are assumed to be in the far field zones of each other and will thus scatter the light independently.
For a snow layer, the snow grains can be located in each other's near-field, resulting in interactions between the scattered electromagnetic fields of neighboring particles which leads to modification of single-scattering properties (Mishchenko, 2014;Mishchenko, 1994). The impact of near-field effect was investigated in Pohl et al. (2019) using the modified single scattering properties of sparsely distributed particles as suggested in Mishchenko (1994). The comparison of snow BRFs calculated assuming sparsely or densely packed snow layers shows that the maximum difference does not exceed 0.015% (Pohl et al., 2019). Therefore, this effect was ignored in radiative transfer calculations through the snow layer.
To account for atmospheric effects, SCIATRAN incorporates a comprehensive database containing monthly and zonal vertical distribution of trace gases e.g. O3, NO2, SO2, H2O, etc., spectral characteristics of gaseous absorbers, vertical profiles of pressure and temperature and molecular scattering characteristics (see Rozanov et al., (2014)  scattering and absorption by aerosols over snow in SCIATRAN, the optical characteristics of aerosol particles and vertical distribution of aerosol number density are required. In this study we use Moderate Resolution Imaging Spectrometer (MODIS) collection 5 aerosol parameterization (Levy et al., 2007) as an internal database in SCIATRAN. Levy et al. (2007) developed a framework for connecting the aerosol micro-physical properties such as the refractive index and size distribution to the AOT at 0.55μm. Using AOT from ground based measurements of AERONET at Barrow/Utqiaġvik as mentioned and selecting one of aerosol types, the Mie code incorporated into SCIATRAN is employed to calculate aerosol extinction and scattering coefficients. In this study, the vertical profile of aerosol number density as an "exponential vertical distribution" for a height of 3.0 km is used.

Sensitivity of reflectance to the snow morphology and atmospheric parameters
The measured reflectance in the visible and NIR spectral range over a snow field depends on the relative importance of the absorption and scattering radiative transfer processes in the atmosphere and snow layer. In this section, we investigate the sensitivity of the reflectance on the radiative transfer through the atmosphere and the snow at the selected wavelength bands: i) 1.649 μm because of the high sensitivity of this wavelength to snow grain properties; and ii) 0.677 and 0.873 μm wavelengths because of the relatively high and differing sensitivities at these wavelengths to the atmospheric conditions and being used for aerosol optical thickness retrievals.

Impact of snow: size and shape of ice crystals
To study the influence of ice crystal morphology on the radiation field above snow covered surfaces, we perform the simulation at 1.649 μm for three important reasons (Leroux et al., 1998): i) the absorption of ice crystals is small or negligible at the selected wavelengths in the visible domain of the spectrum. In contrast, in the near-infrared range, due to the large absorption of ice crystals at these wavelengths, the snow reflectance is significantly affected by the snow grain size; the larger the particle, the smaller the reflectance because of larger absorption and stronger forward scattering; ii) the BRF properties of snow at 1.649 μm are closer to that for single scattering behavior and it is linked to the phase matrix which strongly depends on the shape of ice crystals; iii) the impact of the atmosphere (absorption by CO2 and H2O and diffuse incident irradiance) at 1.649 μm is negligible.
To illustrate the high sensitivity of radiation field to the varying size of ice crystals at 1.649 μm, we simulated snow reflectance at principal and cross planes assuming nine ice crystal morphologies with varying sizes (here size refers to maximum dimension/edge length) 60 ~ 10000 μm and three different roughness (smooth surface: 0, moderate surface roughness: 0.03, severe surface roughness 0.5), for further information see Yang et al. (2013). Fig. 3 shows the simulated reflectance versus the VZA in the principal plane (as the most sensitive and representative direction for the largest changes of reflectance) using severely roughened morphology. As can be seen in Fig. 3, the reflectance strongly changes with the size of ice crystals from 60 to 10000 μm. Differentiating between various shapes has the largest effect in forward scattering (φ=0°) and lesser effect in backward scattering direction (φ=180°). The results indicate that the effect of changing size is larger than the impact of differentiating between various shapes of ice crystals at this wavelength.
Using the "aggregate of 8 columns" shape and changing maximum dimension from 60 μm to 10000 μm result in reflectance decrease of ~ 40 % at nadir (VZA ~ 0°) and more than 80 % in forward scattering direction (at VZA of 60°) which is considerably large. Changing the shape to the droxtal at the same size, provides a noticeable decrease of ~ 30% in reflectance at forward scattering direction for a viewing zenith angle of 60° and leads to a much weaker forward peak.
Noteworthy is, that the plate shape cannot reproduce the enhancement in backward direction (typical for a BRF over snow).
Using aggregate of 5 and 10 plates leads to larger reflectance in all directions. However, the analysis of simulation results at cross plane (not shown here) indicates that, the impact on the reflectance pattern, originating from the specific shapes of the ice crystals is relatively small compared to the impacts at principal plane.
The large range of changes of the reflectance when using different ice crystal sizes in both the principal and cross planes highlights the importance of having accurate priori knowledge or estimation of size of the ice crystals and their shapes to simulate accurately measurements. In our study, due to the lack of such information from in situ measurements, we estimate the size of ice crystals for each selected crystal shape separately to have a priori knowledge of ice crystal properties and limit the differences between the simulated and measured reflectance. The detailed explanation and results are given in sect. 6.

Impact of atmosphere: scattering and absorption by aerosol and gases
The incident radiation on the snow layer is composed of direct sunlight and the diffuse radiation from the sky (Aoki et al., 1999). To take the atmospheric absorption and scattering into account, we assume an atmosphere over the snow layer, which contains: i) Rayleigh scattering (scattering by air molecules), ii) gaseous absorption and iii) absorption and scattering by aerosols. Therefore, in this section, absorption bands e.g. 0.677 μm are selected to evaluate the impact of atmosphere. We calculate the reflectance at 0.677 μm under three different conditions, assuming a model atmosphere governed: i) by Rayleigh scattering; ii) identical to i) but with absorption by ozone (O3) and nitrogen dioxide (NO2); iii) identical to ii) but including aerosol. The calculations are performed assuming following properties of the atmosphere and snow layer: i) Vertical profile of nitrogen dioxide, pressure and temperature are selected according to a 2D chemical transport model (Sinnhuber et al., 2009)

incorporated in SCIATRAN;
ii) Total vertical column of ozone as well as AOT are set according to Table 1; iii) Snow layer is composed of ice crystals having the shape "aggregate of 8 column", maximum dimension of 650 μm and severely roughened crystal surface. gaseous absorption is the smallest ~ 5% close to the nadir region and becomes larger ~ 10% in forward scattering direction which decreases to ~ 8% at 1700 m altitude. At this wavelength, ozone with vertical optical depth (VOD) of 1.6×10 -2 has a much larger contribution to gaseous absorption as compared to that of NO2 with VOD of 3.95×10 -5 .
The reflectance for an atmosphere containing three types of aerosol (weakly/moderately/strongly absorbing aerosol) and without aerosol (containing only Rayleigh and gaseous absorption) are presented in Fig. 4. For more information on aerosol typing used in this study please see Levy et al. (2007). The changes in reflectance due to weakly absorbing aerosol with an AOT of 0.11 (measured by AERONET) at 206 m flight altitude are ~ 5% at nadir and increase in forward scattering direction to ~ 13%. The strongly absorbing aerosol (at the same AOT of 0.11) reduces the reflectance by ~7% at nadir and 2 0 % in forward scattering direction. At 1700 m the reflectance decreases by 6% at nadir and 7% in forward scattering direction. The differences between the three aerosol types does not lead to changes in reflectance, which are larger than 5% in or close to nadir areas. In summary, an atmosphere containing Rayleigh scattering, absorption by ozone (O3) and nitrogen dioxide (NO2) and weakly absorbing aerosol is the best representation of the atmospheric conditions for our case study.

Retrieval of snow grain size and shape
In the previous section, we showed that having adequate a priori information about snow surface and atmosphere is necessary to calculate reflectance of sufficient accuracy. In contrast to the atmospheric parameters available from independent sensors and models, a priori knowledge about ice crystal size and shape for the underlying snow layer is not typically available. To estimate the optimal ice crystal morphology we used a snow grain size and shape retrieval algorithm, by minimizing the difference between the measured and simulated reflectance (See appendix A for details). Here size refers to effective radius 1 of the ice crystal. The retrieval algorithm is applied to measurements at principal and cross planes at 1.649 μm assuming different shape and crystal surface roughnesses. To find the best representative shape and size, the bias and Root Mean Square Error (RMSE) between the measured and simulated reflectance were determined for each case study. Based on comparison, one can state that the angular reflectance pattern of the CAR measurement on the 7 th of April 2008 at Elson Lagoon is reproduced by SCIATRAN successfully. The highest accuracy is obtained by assuming ice crystals as "an aggregate of 8 columns" with severely roughened crystal surface at an effective radius of 98.8 μm (corresponding to maximum dimension of 650 μm). In this case, the largest and smallest discrepancies appear in the forward scattering direction and close to nadir (VZA< ± 25°), respectively. The overall RMSE and bias between measurements and simulation at principal plane is 6.9 % and 2.7 % respectively. A lesser degree of agreement between simulated results and measurements are provided by using "column" and fractal shapes with an RMSE of 7.3% and 9.75%, respectively. The largest difference between measurements and simulations is observed for the case using "droxtal" shape with an RMSE 2 5.54%.
We also retrieved effective radius of ice crystal using CAR data for fresh fallen snow on the 15 th of April 2008. Due to the existing surface horizontal inhomogeneity for the case of fresh snow acquired at lower flight altitude ~ 181 m, larger differences between simulated and measured reflectance are expected, as compared to the old snow case on the 7 th of April 2008. The results are shown in Fig. 6. Unlike the old snow case presented in Fig. 5, the "aggregate of 8 columns" shape does not optimally represent the ice crystals of this particular day. Rather, a reflectance simulated by using an "aggregate of 5 plates" as the ice crystal shape provides the minimum RMSE ~ 12.85% between measurement and simulation results.
"Aggregate of 10 plates" and fractal provide the second and third most representative shapes with an RMSE of ~ 13.16 % and 14.69 %, respectively. The results obtained by using the "droxtal" ice crystal shape exhibit large differences in both of forward and backward scattering directions with RMSE of 34.1 %.
Though the real nature of ice crystal shape at the time of measurement is not known to us, the impact of temperature on morphology of snow grain particles has been debated in previous studies. The results of such studies are now compared with our findings (Slater and Michaelides, 2019;Shultz, 2018;Libbrecht, 2007;Bailey and Hallett, 2004;Yang et al., 2003).
Based on the relationship between temperature and snow grain morphology, the column-based shapes are the dominant ice crystal morphology in environments with temperatures higher than -10°C whereas plates are dominant if the temperature is less than -10°C. Our findings with respect to the most representative shape for each case study agree with this argument. The temperature range during CAR measurements at 6-7 th of April 2008 is from -20 to -5°C. Based on our results the "aggregate of 8 columns" is the most representative shape for measurements conducted on this day. On 15 th of April 2008 when the temperature range changes to -23 to -17°C, mainly plate-based ice crystal shapes are expected for such low temperatures and our results confirm this argument. In addition, the existence of droxtal ice crystals during the measurement is less probable because very low temperatures (~ -50°C) are needed to form droxtal or quasi-spherical ice crystals (Yang et al., 2003). The temperature dependence of the ice crystal morphologies explains in part why droxtal shaped ice crystals do not capture the derived snow reflectance values from CAR measurements in any of our scenarios. With respect to size of ice crystals, we do not compare fresh and old snow cases because it is important to note that the date of old snow case is before fresh snow. This means the studied old snow case is not the aged fresh snow case. Therefore, the change of ice crystal size with its age is not studied in the scope of this paper.
A summary of retrieved effective radii using different ice crystal shapes and corresponding bias and RMSE values is presented in Table 2. The ice crystals with minimum RMSE value at 1.649 μm are underlined and selected to be used for subsequent calculations of reflectance at 0.677, 0.873 μm. In Fig. 7, the importance of ice crystal shape selection for the snow grain size retrieval and the snow reflectance calculation is highlighted. The measurements were selected from the old snow and fresh snow cases. The effective radius is retrieved only at 1.649 μm and then has been used to calculate the reflectance at 0.677 μm and 0.873 μm. The results are presented in Fig. 7 with corresponding RMSE and bias values in the principal plane. It can be seen that the retrieved effective radius value changes from shape to shape. The difference in retrieved effective radius generally does not exceeds 40 % but in the case of the plate ice crystals the retrieved effective radius is ~70 % smaller than the other shapes e.g. aggregate of 8 columns. This is a significantly large difference. However, these results are presented for the principal plane where the maximum differences between simulation and measurement is expected. Therefore, the overall bias and RMSE value on all azimuth direction is smaller than presented here. It can be seen that the RMSE values at 0.677 μm and 0.873 μm are significantly smaller than that at 1.649 μm. This is explained by the high reflectance values at these wavelengths and therefore the larger denominator in RMSE formula, in which the difference of measured and simulated reflectance is divided by measured reflectance.

Comparison of measured and simulated reflectance
In this section we present results of the comparison of measured and simulated reflectance in the snow-atmosphere system.
The simulations, which used the results and findings described in the based on the previous section were performed: assuming an atmosphere containing O3, NO2, weakly absorbing aerosol as described in Table 1. The snow layer is assumed to be comprised of "aggregate of 8 column" ice crystals with a maximum dimension of 650 μm (effective radius 98.83 μm) for the case of old snow, and "aggregate of 5 plates" ice crystals with a maximum dimension of 725 μm (effective radius 83.41 μm) for the case of fresh snow. To assess the accuracy of simulations over all azimuth angles, the correlation plot and the Pearson correlation coefficient between measured and modeled reflectance are shown in Fig. 10.
In Fig. 8, the difference between the simulated and measured reflectance at 0.677 and 1.649 μm is small on average, being less than 0.025 in regions of small VZA and not exceeding ±0.05 for larger VZA < 50°. These values are larger at 0.873 μm; the maximum difference reaches ~ ±0.05 for small VZA. The difference between SCIATRAN simulation values and those of the measurements is pronounced in the forward scattering region where |∆φ| < 40°. As it is shown in Fig. 10, the correlation coefficient between reflectance measurements over old snow and simulation is high, ~ 0.98. Fig. 9 is the same plot as Fig. 8 but for fresh snow. We consider that surface inhomogentities and related larger shadowing effects at measurement altitude of 181 m, explain why the correlation decreased to 0.88 at 0.677 μm (for the case of old snow, acquired at 647 m flight altitude, surface inhomogenities are smoothed and therefore the old snow case is less affected by surface inhomogenities). However, in Fig. 10, at 1.649 μm correlation coefficient is high ~ 0.97 for the case of fresh snow, possibly because of their being less sensitivity of this channel to shadowing and atmospheric scattering effects. In addition, the high correlation coefficient at 1.649 μm and small discrepancies < ±0.025 in off-glint region confirms the suitability of the selection of the best representation for ice crystal shape in previous step. The differences between SCIATRAN simulations and CAR measurements of reflectance are less pronounced in the glint region, as compared to those for the old snow.

Conclusion
In this study, our objective was to assess the accuracy of the simulation of the reflectance in a snow-atmosphere system taking different snow morphology and atmospheric absorption and scattering into account. For this we used a state of the art RTM, SCIATRAN, which used explicit models of the snow layer and the airborne CAR measurements. The SCIATRAN RTM (a phenomenological RTM) was used to simulate the reflectance in the snow-atmosphere system and its changes for different snow morphologies (i.e. snow grain size and shape). These simulations take atmospheric scattering and absorption explicitly into account. We investigated the sensitivity of reflectance in the snow-atmosphere system to snow grain size and shape. We have shown that the selection of the most representative shape and size of the nine ice crystals used in SCIATRAN to describe the snow surface is essential to minimize the difference between simulations and measurements.
To obtain a priori knowledge of snow morphology, we use the snow grain size and shape retrieval algorithm and apply it to CAR data. In our case study at Barrow/Utqiaġvik, the simulated reflectance assuming ice crystals with aggregate composed of 8 columns shape agreed well with measurements for the old snow case, having RMSE of 6.9 % and average bias of 2.7 % with respect to the measured CAR reflectance in the principal plane where the largest discrepancies are expected. For the case of freshly fallen snow, an aggregate of 5 plates shape was the most representative ice crystal having RMSE values of 12.8 % and a bias of 11.23 % with respect to the measured CAR reflectance. The data for the freshly fallen snow case were acquired at 181 m. Larger differences as compared to the older snow case at 647 m are attributed to surface inhomogenity. The surface inhomogenity most likely originate from sastrugi. Simulation, in which the snow layer was comprised of ice crystals with a droxtal shape (being semi-spherical particles) did not yield accurate reflectance for the snow-atmosphere system in any of our case studies. We showed that using the knowledge from studies of the temperature dependence of ice crystal morphologies agrees with our findings with respect to the most representative ice crystal size and shape for our case studies.
In our study, the simulated patterns of the reflectance with respect to spectral and directional signatures produce well the measurements, as evidenced by the high correlation coefficients in the range of 0.88 ~ 0.98 between measurements (old and fresh snow) and simulation at the selected wavelengths of 0.677, 0.873 and 1.649 μm. In the off-glint regions |∆φ| > 40° and VZA < 50°, the overall absolute difference between the modeled reflectance from SCIATRAN and CAR measurements is below 0.05. This absolute difference in off-glint area is smaller in the short wave infrared as compared to visible. It should be noted here that the reflectance of the snow is lower in the short wave infrared compared to the visible.
In summary, the approach shows the high accuracy of the phenomenological SCIATRAN RTM in simulating the radiation field in the snow-atmosphere scenes for off-glint observations. The results are applicable for the inversion of snow and atmospheric data products from the satellite or airborne passive remote sensing measurements above snow. To mitigate the relatively larger differences between measurements and simulation for glint condition as compared to off-glint region, the use of a vertically inhomogeneous snow layer consisting of different ice crystal shapes and sizes is proposed. This research has been undertaken as part of the investigations in the framework of trans regional (AC)3 project (Wendisch et al., 2017) that aims to identify and quantify different parameters involved in rapidly changing climate in the Arctic. In this respect, the analysis of this study will be used to improve the assumptions made for reflectance in snowatmosphere system in the algorithms designed to retrieve atmospheric parameters (such as AOT) above Polar Regions.

Appendix
For the selected snow models using different ice crystal morphologies, the variation of the snow reflectance R(λ,Ω) at wavelength λ and direction Ω with respect to the variation δre(z) of the effective radius profile re(z) along the vertical coordinate z within snow layer can be presented, neglecting nonlinear terms, by the following equation: where R 0 ( λ , Ω ) and R ( λ , Ω) are the reflection functions calculated assuming an effective radius profile ŕ e ( z ) and ŕ e ( z )+δ r e ( z) , respectively. The angular variable Ω={θ 0 ,θ , φ} comprises a set of variables: θ 0 is the solar zenith angle, θ and φ are the zenith and relative azimuthal angles of observation direction; Z t is the top altitude of snow layer and is the functional derivative of the function R ( λ , Ω) with respect to the function r e ( z ) which is also called weighting function (Rozanov et al., 2007). The weighting function was calculated using a numerically efficient forward-adjoint approach (Rozanov, 2006;Rozanov and Rozanov, 2007) implemented in the SCIATRAN model. Here, it is assumed that properties of snow do not change in the horizontal plane and within snow layer there is no additional absorber such as soot, dust or other pollutants. We note that the weighting function includes the contribution of variations not only by the scattering and extinction coefficients but also by the phase function.
Although the linear relationship given by Eq. (1) can be used to retrieve the vertical profile of the effective radius within the snow layer in a way similar to that used to the morphology of water droplets (Kokhanovsky and Rozanov, 2012), we restrict ourselves to the assumption of independent of the altitude r e . Introducing the weighting function for the absolute variation of the effective radius as: we have The resultant linear relationship is a basic equation to formulate inverse problem with respect to the parameter r e using measurements of spectral reflectance.
For practical applications Eq. (4) should be rewritten in the vector-matrix form as follows: 14 The components of vectors Y and Y 0 are the measured and simulated reflectance at discrete number of observation directions Ω j and wavelengths λ i , the elements of matrix K are weighting functions W r ( λ i , Ω j ), X=[r e ] is the state vector, is the a priori state vector. We note that in the case under consideration, the matrix K and state vector X are represented by the column vector and scalar, respectively.
Assuming that the number of discrete observation directions Ω j and wavelengths λ i is larger than the dimensions of the state vector, the solution of Eq. (5) is obtained by minimizing the following cost function: which describes the root-mean-square deviation between measured and simulated snow reflectance.
Owing to the linear relationship given by Eq. (5) the minimization problem formulated above can be solved analytically: In deriving Eq. (7) we have neglected the linearization error which can be significant if X 0 is far from X. To mitigate the impact of linearization error we solve minimization problem given by Eq. (6) iteratively. In particular instead of Eq. (7) is used X n = X n −1 +(K n −1 T K n −1 ) Where n=1, 2, … is the iteration number, K n −1 and Y n −1 are the matrix of weighting functions and reflectance vector calculated using the state vector X n −1 . The iteration process is finished if the difference between X n and X n −1 is smaller than a preselected criteria.
The calculation of weighting functions and TOA reflectance is performed at each iteration step using the radiative transfer model SCIATRAN (Rozanov et al., 2014). In SCIATRAN weighting functions are calculated employing a very efficient forward-adjoint technique, which is based on the joint solution of the linearized forward and adjoint radiative transfer equations (Rozanov 2006;Rozanov and Rozanov 2007 and references therein). This enables the TOA reflectance and required weighting function to be calculated simultaneously.   The green lines indicate simulated reflectance assuming Rayleigh scattering (case i); the blue line shows reflectance for case ii (as case i including absorption of O3 and NO2, the orange lines show the reflectance for case iii (as case ii but adding aerosol with an AOT of 0.11 for three types of aerosol: i) weakly absorbing, ii) moderately absorbing and iii) strongly absorbing).