Simulated reflectance above snow constrained by airborne measurements of solar radiation: implications for the snow grain morphology in the Arctic

Accurate knowledge of the reflectance from snow/ice-covered surfaces is of fundamental importance for the retrieval of snow parameters and atmospheric constituents from space-based and airborne observations.

In this paper, we simulate the reflectance in a snow–atmosphere system, using the phenomenological radiative transfer model SCIATRAN, and compare the results with that of airborne measurements. To minimize the differences between measurements and simulation, we determine and employ the key atmospheric and surface parameters, such as snow grain morphologies (or habits).

First, we report on a sensitivity study. This addresses the requirement for adequate a priori knowledge about snow models and ancillary information about the atmosphere. For this aim, we use the well-validated phenomenological radiative transfer model, SCIATRAN. Second, we present and apply a 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 for estimating the most representative snow model, using different types of snow morphologies, for the airborne observation conditions performed by NASA's Cloud Absorption Radiometer (CAR).

Third, we present a comprehensive comparison of the simulated reflectance (using retrieved snow grain size and shape and independent atmospheric data) with that from airborne CAR measurements in the visible (0.670 µm) and near infrared (NIR; 0.870 and 1.6 µm) 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 the snow layer consists of ice crystals with aggregates of eight column ice habit and having an effective radius of ∼99 µm, we find that, for a surface covered by old snow, the Pearson correlation coefficient, R, between measurements and simulations is 0.98 (R²∼0.96). For freshly fallen snow, assuming that the snow layer consists of the aggregate of five plates ice habit with an effective radius of ∼83 µm and having surface inhomogeneity, the correlation is ∼0.97 (R²∼0.94) in the infrared and 0.88 (R²∼0.77) in the visible wavelengths. The largest differences between simulated and measured values are observed in the glint area (i.e., in the angular regions of specular and near-specular reflection), with relative azimuth angles $<\pm \mathrm{40}$^∘ in the forward-scattering direction. The absolute difference between the modeled results and measurements in off-glint regions, with a viewing zenith angle of less than 50^∘, is generally small $\sim \pm \mathrm{0.025}$ and does not exceed ±0.05. These results will help to improve the calculation of snow surface reflectance and relevant assumptions in the snow–atmosphere system algorithms (e.g., aerosol optical thickness retrieval algorithms in the polar regions).

1 Introduction

The extent and type of snow and ice cover have a significant impact on climate, as noted by Arrhenius over 100 years ago (Arrhenius, 1896). There is a positive feedback between decreasing surface temperature, an increase in snow and ice cover, and an associated increase in planetary albedo, which then further decreases surface temperature and vice versa. Consequently, changes in snow and ice extent and morphology play a key role in climate change. Having accurate knowledge about snow/ice-covered surfaces is a prerequisite for identifying and quantifying changes in the climate (Schneider and Dickinson, 1974; Curry et al., 1995; Cohen et al., 2014; Kim et al., 2017; Wendisch et al., 2017, 2019).

In addition, during recent decades, the Arctic region has warmed more rapidly than other regions. This phenomenon is known as the Arctic amplification (AA; Serreze and Barry, 2011). The analysis of the growing number of long-term records of the data products (e.g., the amount of trace gases, aerosol, and cloud parameters) retrieved from passive and active satellite observations provides potentially invaluable information for identifying and quantifying the evolution and consequences of AA (Wendisch et al., 2017).

Because of the magnitude of the scattering from snow, the use of remote sensing measurements above snow-covered surfaces in the cryosphere requires accurate models of the scattering and reflectance from snow surfaces. However, the current differences between simulated and measured reflectance in a coupled snow–atmosphere system lead to systematic errors in the determination of the atmospheric constituents, in particular for clouds and aerosol parameters but also trace gases (e.g., Istomina et al., 2010; Istomina, 2012; Jafariserajehlou et al., 2019).

A large number of experimental and theoretical studies have been conducted measuring and modeling snow optical properties, such as the angular distribution of reflected light over and within the snow surface and the subsequent derivation of snow albedo. The early measurements by Middleton and Mungal (1952) and the model of Dunkle and Bevans (1956) used to analyze the transmittance and reflectance of snow cover were the beginning of considerable efforts on this topic. Barkstrom (1972) formulated and solved the scattering problem for snow surfaces in terms of the radiative transfer theory. Later, the comparison of reflectance from snow-covered surfaces simulated by radiative transfer models (RTMs) with that of observation made substantial progress in our understanding of the angular distribution of reflectance above snow surface (e.g., Wiscombe and Warren 1980; Warren et al., 1998; Arnold et al., 2002; Painter et al., 2003; Kokhanovsky and Zege, 2004; Li and Zhou, 2004; Hudson et al., 2006; Hudson and Warren, 2007; Lyapustin et al., 2010; Kokhanovsky and Breon, 2012).

The reflection/scattering patterns of snow-covered surfaces can be summarized as follows: (i) snow is not a Lambertian reflector in the visible and near infrared spectral region because its reflectance has an anisotropic nature, and the anisotropy increases with wavelength; (ii) unlike other surface types (e.g., vegetation or soil) with a strong peak in back-scattering direction (the hot spot effect), snow has a strong forward peak for large viewing zenith angles (e.g., Gatebe and King, 2016); (iii) the snow reflectance variation is larger in the principal plane (i.e., the plane containing the Sun, surface normal, and observation direction) than in the cross plane (i.e., the one perpendicular to principal plane; Warren, 1982; Lyapustin et al., 2010; Kokhanovsky and Breon, 2012). However, the remaining discrepancies between simulated reflectance in a snow–atmosphere system and field measurements led to further investigations in the field of the single scattering properties of snow grains (Mishchenko et al., 1999; Jin et al., 2008; Yang and Liou, 1998; Yang et al., 2003, 2013), surface roughness (Warren et al., 1998; Hudson et al., 2006; Hudson and Warren, 2007; Lyapustin et al., 2010; Zhuravleva and Kokhanovsky, 2011), and atmospheric-correction methods (Lyapustin et al., 2010). Despite substantial improvements, the uncertainties in our understanding of the microphysical and macroscopic properties of snow are an unresolved issue for RTMs, ray tracing, and climate models. For example, the current state-of-the-art RTMs yield much more anisotropic reflectance behavior for snow in the glint region than observed in reality (Zhuravleva and Kokhanovsky, 2011; Lyapustin et al., 2010; Hudson and Warren, 2007; Warren et al., 1998). These studies either focus on the snow reflectance at the surface, employing an atmospheric-correction method (Leroux et al., 1998; Kokhanovsky and Zege, 2004; Kokhanovsky et al., 2005; Lyapustin et al., 2010; Negi and Kokhanovsky, 2011), or consider the atmospheric effects without in-depth investigations of the surface parameters (Aoki et al., 1999; Hudson et al., 2006; Kokhanovsky and Breon, 2012). A comprehensive study and investigation of both the snow layer and atmosphere parameters in a coupled snow–atmosphere system has not yet been undertaken but is required to improve the accuracy of remote sensing retrieval algorithms for aerosol and cloud in the Arctic region (Istomina et al., 2010; Istomina, 2012; Jafariserajehlou et al., 2019).

The intent of this paper is to (i) study the sensitivity of scattering and reflectance in the coupled snow–atmosphere system to both surface and atmospheric parameters, (ii) retrieve the most representative ice crystal morphology by applying a snow grain size and shape retrieval algorithm to measured reflectance, and (iii) evaluate the ability of a phenomenological RTM to reproduce the measured reflectance over the spectral range 0.34–1.649 µm at all available observation directions using the retrieved atmospheric and snow parameters.

In this study, we use the RTM SCIATRAN (Rozanov et al., 2014), which is a well-validated phenomenological RTM, and the airborne observations of reflectance acquired by the Cloud Absorption Radiometer (CAR). The CAR measurements were made during the Arctic Research of the Composition of the Troposphere from Aircraft and Satellite (ARCTAS) spring 2008 campaign over Barrow/Utqiaġvik, Alaska. Further information about the atmospheric parameters during the ARCTAS campaign was taken from available Aerosol Robotic Network (AERONET) and satellite data.

The paper is organized as follows: in the next section, we present the theoretical background and terminology used to calculate the angular distribution of reflectance in a snow–atmosphere system. In Sects. 3 and 4, we introduce and explain the measurements and the simulation methods. In Sect. 5, the sensitivity of reflectance to the underlying snow layer and atmospheric parameters is investigated. In Sect. 6, the results of applying the two-stage snow grain size and shape 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 A contains a detailed description of the snow grain size and shape retrieval algorithm used in the study.

2 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 follows (Nicodemus et al., 1977; Schaepman-Strub et al., 2006):

$$\begin{array}{}\text{(1)}& {\text{BRDF}}_{\mathit{\lambda }}={\displaystyle \frac{\mathrm{d}{L}_{\mathrm{r}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{r}},{\mathit{\phi }}_{\mathrm{r}};\mathit{\lambda }\right)}{\mathrm{d}{E}_{\mathrm{i}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}};\mathit{\lambda }\right)}}\left[{\mathrm{sr}}^{-\mathrm{1}}\right],\end{array}$$

where L_r is the reflected radiance, and θ and φ are the zenith and azimuth angles, respectively. The subscripts of i and r correspond to the incident and to the reflected beams, respectively. E is the incident surface flux (irradiance), and λ is the wavelength. However, the BRDF is not a directly measurable quantity because of its formulation 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 follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{e}}& ={\displaystyle \frac{\mathrm{\Delta }{L}_{\mathrm{r}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{r}},{\mathit{\phi }}_{\mathrm{r}};\mathit{\lambda }\right)}{\mathrm{\Delta }{E}_{\mathrm{i}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}};\mathit{\lambda }\right)}}\\ & ={\displaystyle \frac{\mathrm{\Delta }{L}_{\mathrm{r}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{r}},{\mathit{\phi }}_{\mathrm{r}};\mathit{\lambda }\right)}{\mathrm{\Delta }{L}_{\mathrm{i}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\phi }}_{\mathrm{i}};\mathit{\lambda }\right)\cos {\mathit{\theta }}_{\mathrm{i}}\mathrm{\Delta }{\mathit{\omega }}_{\mathrm{i}}}}\left[{\mathrm{sr}}^{-\mathrm{1}}\right],\end{array}\end{array}$$

where ${\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{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, subsurface scattering, and radiation parameters such as wavelength, polarization, and fluorescence (i.e., significant variations do not occur within small intervals; see Nicodemus et al., 1977; Gatebe and King, 2016). As a result, the ${\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{e}}$ is determined by the following:

$$\begin{array}{}\text{(3)}& {\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{e}}={\displaystyle \frac{L_{\mathrm{r}}^{\mathrm{e}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{r}},\mathrm{\Delta }\mathit{\phi }\right)}{F_{\mathrm{0},\mathit{\lambda }}\cos {\mathit{\theta }}_{\mathrm{i}}}},\end{array}$$

where $L_{\mathrm{r}}^{\mathrm{e}}$ is the measured radiance, and 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, it is the equivalent bidirectional reflectance factor ${\text{BRF}}_{\mathit{\lambda }}^{\mathrm{e}}$, for which ${\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{e}}$ multiplied by π is more representative.

To isolate the reflectance properties of the surface and derive ${\text{BRF}}_{\mathit{\lambda }}^{\mathrm{e}}$ or ${\text{BRDF}}_{\mathit{\lambda }}^{\mathrm{e}}$ just above the surface, we need to apply atmospheric-correction methods on the measured radiance at TOA or flight altitude (e.g., by using knowledge of the atmospheric scattering or absorption by applying RTMs). This removes the following four atmospheric contributions from the measured radiance at TOA or flight altitude (Schaepman-Strub et al., 2006): the contribution of light scattered by the atmosphere (i) before the solar radiation has reached the surface, (ii) after being reflected by the surface, (iii) before and after reaching the surface, and (iv) the atmospheric path radiance.

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 the wavelength, in the range of 0.381 to 2.324 µm (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).

The atmospheric-correction method relies on different assumptions by which several sources of uncertainty should be taken into account. In this study, to avoid such uncertainties, we do not apply an atmospheric correction to the measurements (radiances L_r, h) at flight altitude (h). Instead, we calculate and use the reflectance at flight altitude with the following equation:

$$\begin{array}{}\text{(4)}& R={\displaystyle \frac{\mathit{\pi }{L}_{\mathrm{r},\mathrm{h}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{r}},\mathrm{\Delta }\mathit{\phi }\right)}{F_{\mathrm{0},\mathit{\lambda }}\cos {\mathit{\theta }}_{\mathrm{i}}}},\end{array}$$

where L_r, h is the measured radiance at flight altitude. All reflectance/${\text{BRF}}_{\mathit{\lambda }}^{\mathrm{e}}$ values at flight altitude in this study represent R in Eq. (4) and are referred to as the reflectance factor in the snow–atmosphere system.

In the simulation of the reflectance factor in a coupled snow–atmosphere system, we need to account for the atmospheric effects contribution properly. For this reason, we take independent data about atmospheric parameters (aerosol optical thickness, AOT, and gases absorption) from ground-based and space-borne measurements. We select the data with the closest spatial and temporal interval actual airborne measurements. We discuss more details of the atmospheric data and their application to the simulation routine in Sects. 3 and 4. To estimate ${\text{BRF}}_{\mathit{\lambda }}^{\mathrm{e}}$ just above the surface, further atmospheric correction is needed. We assume the reflectance factor at flight altitude is a good estimation of ${\text{BRF}}_{\mathit{\lambda }}^{\mathrm{e}}$ just above the surface at infrared wavelengths where atmospheric scattering is negligible.

3 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 from 1984 up to the present. CAR has been used to measure the single scattering albedo of clouds, the bidirectional reflectance of various surface types, and for acquiring imagery of clouds and the Earth's surface. For this study, we used CAR data from the ARCTAS campaign conducted at Elson Lagoon, near Barrow/Utqiaġvik, Alaska, in April 2008 as part of the International Polar Year (Lyapustin et al., 2010; Gatebe and King, 2016). The goal of ARCTAS was to study physical and chemical processes in the Arctic atmosphere and (e.g., long-range transport of pollution to the Arctic) and surface parameters (e.g., snow reflectance angular variation). A P-3B aircraft carried the CAR instrument and was deployed by NASA from Fairbanks, Alaska. Figure 1 shows the flight track on 7 April 2008. The date, location, measurement geometry, and available atmospheric parameters during the measurements used in this study are presented in Table 1.

Figure 1 Flight track of P-3B airplane carrying NASA's Cloud Absorption Radiometer (CAR) on 7 April 2008 during the Arctic Research of the Composition of the Troposphere from Aircraft and Satellite (ARCTAS) campaign (figure adapted from: https://www-air.larc.nasa.gov/missions/arctas/arctas.html, last access: November 2020). Image credits: © NASA; © Tele Atlas 2008; © Europa Technologies 2008; © TerraMetrics 2008; © 2007 Google.

Table 1 Summary of the CAR, Aerosol Robotic Network (AERONET) aerosol optical thickness (transferred from 0.5 to 0.55 µm), and weighting function differential optical absorption spectroscopy (WFDOAS) ozone data used in this study. Note: solar zenith angle – SZA.

The unique design of CAR simultaneously provides both up-welling and down-welling radiances at 14 spectral bands (Table 2) located in the atmospheric window regions of UV, visible, and near-infrared from 0.34 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 the data with a mirror rotating 360^∘ in a plane perpendicular to the direction of flight through a 190^∘ aperture that allows the acquisition of 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 the viewing zenith and azimuth angles allows the estimation of the anisotropy of the reflectance in the snow–atmosphere system with high accuracy.

Table 2 Summary of CAR wavelengths and bandwidth.

The spatial resolution of CAR depends on the flight altitude, e.g., 10 and 18 m² at nadir for 600 and 1000 m flight altitude, respectively, which increases with the viewing zenith angle (VZA), e.g., 580 m² at 80^∘ VZA for 1000 m flight altitude. The capability of acquiring data at different altitudes (∼200, 600, and 1700 m) enables us to evaluate the sensitivity of reflectance with respect to atmospheric effects in RTM simulations.

Examples of calculated reflectance factor values using Eq. (4) from CAR measurements on 7 April 2008 at Elson Lagoon (71.3^∘ N, 156.4^∘ W) are shown in Figs. 2 and 3. As we can see in these two figures, in spite of the influence of the atmospheric scattering and absorption, the general features of the snow BRF are clearly observable in polar plots and principal and cross-plane plots because of (i) the decrease in snow reflectance with increasing wavelength due to the increasing absorption by snow at longer wavelengths, (ii) the increase in the snow BRF as a function of VZA and the strong forward-scattering peak in the principal plane at large VZA, and (iii) the smaller angular variation in the BRF at cross plane compared to the principal plane, though the reflectance values increase with VZA. The snow surface spatial inhomogeneity decreases with increasing altitude due to the change in spatial resolution with altitude (Gatebe and King, 2016; Lyapustin et al., 2010). Accordingly, at poorer spatial resolutions, spatial homogeneity is more efficiently averaged, as can be seen in Fig. 2 at flight altitude of 1700 m, compared to 206 m, in which we have a higher spatial resolution.

Figure 2 Angular distribution of reflectance factor in the snow–atmosphere system derived from CAR measurements on 7 April 2008 at Elson Lagoon (71.3^∘ N, 156.4^∘ W), showing reflectance at the 0.677 µm wavelength and three flight altitudes of 206, 647, and 1700 m, respectively (a, b, c), and the 1.032 µm wavelength and the same flight altitudes (d, e, f). The principal plane is the horizontal line (φ=0 and 180^∘), and the viewing zenith angle is shown as the radius of the polar plots from 0^∘ (nadir) to 70^∘. The solar zenith angle is 70.23, 69.11, and 67.78^∘ for flight altitudes of 206, 647, and 1700 m, respectively.

Figure 3 Angular distribution of the reflectance factor in the snow–atmosphere system, derived from measurements by CAR at 647 m flight altitude and six wavelengths of 0.677, 0.873, 1.032, 1.222, 1.275, and 1.649 µm on 7 April 2008 at Elson Lagoon (71.3^∘ N, 156.4^∘ W). (a) The principal plane (φ=0 and 180^∘) and (b) cross plane (φ=90 and 270^∘).

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 sunlight. The data are provided often at a high temporal resolution of 15 min. AERONET AOT data are provided at 0.5 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. The aerosol condition and its chemical and optical properties have been measured continuously at Barrow, Alaska, during different seasonal periods (Quinn et al., 2002). Previous studies indicate the largest contribution from sea salt, non-sea-salt sulfate, and mineral dust. The average contribution of black carbon is very small compared to other aerosol types (Udisti et al., 2020). During the haze season (January to April), sea salt plays the dominant role in controlling light scattering in wintertime and non-sea salt sulfate in spring (Quinn et al., 2002). The increase in non-sea-salt (nss) sulfate in January to May is the long-range transport of anthropogenic primary nss sulfate, besides the long-range transport of anthropogenic SO₂ and its photo-oxidization to nss sulfate with the increase of light levels, and the local production of biogenic nss sulfate.

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 differential optical absorption spectroscopy (WFDOAS) algorithm version 4 (Weber et al., 2018). This data set (covering 1995 to the present) consists of merged total ozone column data retrieved by WFDOAS from the Global Ozone Monitoring Experiment (GOME), Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY), and GOME-2A. In this paper, the ozone total column data are selected using the criteria of having the 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 the 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 snow–atmosphere system.

4 Simulations

SCIATRAN is a software package for radiative transfer modeling developed at the Institute of Environmental Physics (IUP), University of Bremen (Rozanov et al., 2002, 2014), and freely available at http://www.iup.uni-bremen.de/sciatran/ (last access: February 2020). The SCIATRAN package has been used in a variety of remote sensing studies to simulate radiative transfer processes in a wide spectral range, from the ultraviolet to the thermal infrared (0.18–40 µm), assuming either a plane parallel or a spherical atmosphere (Rozanov et al., 2014).

In this paper, to calculate the reflectance factor 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., black carbon) 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 nine ice crystal shapes/habits developed by Yang et al. (2013) has been incorporated in the SCIATRAN model (Pohl et al., 2020). This database comprises a full set of single scattering properties at wavelengths from the UV to the far infrared (IR) for the following nine ice crystal morphologies: droxtal, column and hollow column, aggregate of eight columns, plates, small aggregate of five plates, large aggregate of 10 plates, and hollow and solid bullet rosettes. More detailed information about the ice crystal shapes and sizes can be found in Yang et al. (2013). In addition to the abovementioned nine ice crystals, optical parameters for triadic Koch fractal (referred to 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, we use the second-generation fractals, as described in Macke et al. (1996) and Rozanov et al. (2014).

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 lead to the modification of single scattering properties (Mishchenko, 2014; Mishchenko, 1994). The impact of the near-field effect was investigated in Pohl et al. (2020), 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., 2020). 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 distributions of trace gases, e.g., O₃, NO₂, SO₂, H₂O, etc., spectral characteristics of gaseous absorbers, vertical profiles of pressure and temperature, and molecular scattering characteristics (see Rozanov et al., 2014, for details). To account for 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 microphysical 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 the 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.

For the conditions described above, the radiative transfer calculations are performed at a source–target–sensor geometry extracted from the airborne measurements at solar zenith angles of 70.23, 69.11, 67.68, and 62.11^∘ together with a viewing zenith angle (0–70^∘) and relative azimuth angle (0–360^∘), with an angular resolution of 5^∘, and at four different altitudes of 181, 206, 647, and 1700 m. More detailed information about atmospheric and snow layer parameters are given and discussed separately in the following section.

5 Sensitivity of reflectance to the snow morphology and atmospheric parameters

The measured reflectance in the visible and near infrared (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 being used for aerosol optical thickness retrievals.

Figure 4 The change in reflectance factor values in the principal plane (φ=0 and 180) with the size and shape of ice crystals at the wavelength of 1.649 µm. The left column of the legend in each panel shows the maximum length of the ice crystal, and the right column is its equivalent effective radius.

5.1 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 the following 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 those for single scattering behavior, and this is linked to the phase function, which strongly depends on the shape of ice crystals.

• iii.

  The impact of the atmosphere (absorption by CO₂ and H₂O and diffuse incident irradiance) at 1.649 µm is negligible.

To illustrate the high sensitivity of the radiation field to the varying sizes of ice crystals at 1.649 µm, we simulated the snow reflectance factor at principal and cross planes, assuming nine ice crystal morphologies with varying sizes (here size refers to maximum dimension/edge length), 60∼10 000 µm, and three different roughnesses (smooth surface – 0; moderate surface roughness – 0.03; severe surface roughness – 0.50). For further information, see Yang et al. (2013). Figure 4 shows the simulated reflectance factor 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. 4, the reflectance factor strongly changes with the size of ice crystals from 60 to 10 000 µm. The equivalent effective radius¹ is shown beside the maximum lengths of ice crystals in Fig. 4. Differentiating between various shapes has the largest effect in the forward-scattering (φ=0^∘) and a 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 eight columns shape and changing the maximum dimension from 60 to 10 000 µm results in a reflectance decrease of ∼40 % at the nadir (VZA∼0^∘) and more than 80 % in the forward-scattering direction (at VZA of 60^∘), which is considerably large. If we change only the shape of the snow grain from an aggregate of eight columns to the droxtal, but we keep the size (largest dimension) as it is (e.g., 300 µm), this change provides a noticeable decrease of ∼30 % in reflectance at the forward-scattering direction for a viewing zenith angle of 60^∘, and this leads to a much weaker forward peak. What is noteworthy is that the plate shape cannot reproduce the enhancement in a backward direction (typical for a BRF over snow) that is as strong as what the aggregate of eight columns or the droxtal shape causes. Using the aggregate of five and 10 plates leads to larger reflectance in all directions compared to the single plate shape. However, the analysis of the simulation results at the 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 the principal plane.

The large range of changes in the reflectance in both the principal and cross planes, when using different ice crystal morphologies, highlights the importance of having accurate a priori knowledge or estimations of the size and shapes of the ice crystals to accurately reproduce 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 to limit the differences between the simulated and measured reflectance. The detailed explanation and results are given in Sect. 6.

5.2 Impact of atmosphere: scattering and absorption by aerosol and gases

The incident radiation 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 the atmosphere. We calculate the reflectance factor 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 (O₃) and nitrogen dioxide (NO₂), and (iii) identical to (ii) but including aerosol. The calculations are performed assuming the 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 and AOT are set according to Table 1.

• iii.

  Snow layer is composed of ice crystals having the shape aggregate of eight columns, a maximum dimension of 650 µm, and a severely roughened crystal surface.

Figure 5 shows the impact of the atmosphere and the difference between measured and simulated reflectance factor values at three different altitudes (206, 647, and 1700 m) for the three scenarios. The reflectance reduction at 647 m flight altitude due to gaseous absorption is the smallest (∼5 %) close to the nadir region and becomes larger (∼10 %) in the forward-scattering direction, which decreases (to ∼8 %) at 1700 m altitude. At this wavelength, ozone with a vertical optical depth (VOD) of $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{2}}$ has a much larger contribution to gaseous absorption compared to that of NO₂ with a VOD of $\mathrm{3.95}\times {\mathrm{10}}^{-\mathrm{5}}$.

Figure 5 Measured and simulated reflectance factor at 0.677 µm versus the viewing zenith angle (VZA) in the principal plane (φ=0 and 180^∘) at three different flight altitudes. Panels represent results at 206, 647, and 1700 m flight altitudes, respectively. The green lines indicate simulated reflectance assuming Rayleigh scattering (case i). The blue line shows reflectance for case ii (as in case i, including absorption of O₃ and NO₂). The orange lines show the reflectance for case iii (as in case ii but adding aerosol with an aerosol optical thickness (AOT) of 0.11 for three types of aerosol, namely (i) weakly absorbing, (ii) moderately absorbing, and (iii) strongly absorbing).

The reflectance for an atmosphere containing three types of aerosol (weakly, moderately, or strongly absorbing aerosol) and without aerosol (containing only Rayleigh and gaseous absorption) is presented in Fig. 5. For more information on the aerosol typing used in this study 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 the nadir and an increase in the forward-scattering direction to ∼13 %. The strongly absorbing aerosol (at the same AOT of 0.11) reduces the reflectance by ∼7 % at the nadir and ∼20 % in the forward-scattering direction. At 1700 m, the reflectance decreases by 6 % at the nadir and 7 % in the 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 (O₃) and nitrogen dioxide (NO₂), and weakly absorbing aerosol is the best representation of the atmospheric conditions for our case study.

6 Retrieval of snow grain size and shape

In the previous section, we showed that having adequate a priori information about the snow surface and atmosphere is necessary for calculating the reflectance factor with 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 factor (see Appendix A for details). Here, size refers to the effective radius 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 factor were determined for each case study.

Figure 6 shows one example of the comparison between the measured and simulated reflectance factor at principal and cross planes. The absolute uncertainty of CAR measurements is within 5 % and is shown by the uncertainty envelope (shading). The accuracy of our radiative transfer calculations is estimated to be in the range of 0.1 %. Based on the comparison, one can state that the angular reflectance pattern of the CAR measurement on 7 April 2008 at Elson Lagoon is reproduced by SCIATRAN successfully. The highest accuracy is obtained by assuming ice crystals as an aggregate of eight columns with severely roughened crystal surface at an effective radius of 98.8 µm (corresponding to a maximum dimension of 650 µm). In this case, the largest and smallest discrepancies appear in the forward-scattering direction and close to the nadir ($\text{VZA}<\pm \mathrm{25}$^∘), respectively. The overall RMSE and bias between measurements and simulation at the principal plane is 6.9 % and 2.7 %, respectively. A lesser degree of agreement between simulated results and measurements is 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 the droxtal shape with an RMSE ∼25.54 %.

Figure 6 Comparison of measured and simulated reflectance factors. Measurements were performed by the CAR instrument over old snow at 647 m flight altitude on 7 April 2008 at 1.649 µm. The uncertainty in CAR measurements is indicated by the envelope (shading). SCIATRAN simulations in the principal and cross plane given by the dashed–dotted and dotted lines, respectively, and by different colors (green, blue, and red) present smooth, moderately roughened, and severely roughened crystal surfaces. Positive and negative VZAs correspond to azimuthal angles φ=0 and 180^∘ for the principal plane and φ=90 and 270^∘ for the perpendicular plane, respectively.

We also retrieved the effective radius of ice crystal using CAR data for fresh fallen snow on 15 April 2008. Due to the existing surface horizontal inhomogeneity for the case of fresh snow acquired at a lower flight altitude ∼181 m, larger differences between simulated and measured reflectance are expected, compared to the old snow case on 7 April 2008. The results are shown in Fig. 7. Unlike the old snow case presented in Fig. 6, the aggregate of eight columns shape does not optimally represent the ice crystals of this particular day. Rather, a reflectance simulated by using an aggregate of five plates as the ice crystal shape provides the minimum RMSE ∼12.85 % between the measurement and simulation results. An 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 an RMSE of 34.1 %.

Figure 7 The same as Fig. 6 but the measurements by the CAR instrument were performed on 15 April at 181 m flight altitude over fresh snow.

Though the real nature of the ice crystal shape at the time of measurement is not known to us, the impact of temperature and supersaturation on the morphology of snow grain particles has been debated in previous studies (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. Though more investigation is needed, especially to account for the temperature profile at the exact time of snowfall, our findings with respect to the most representative shape for each case study agree with this argument. The temperature range during CAR measurements on 6–7 April 2008 is from −20 to −5 ^∘C. Based on our results, the aggregate of eight columns is the most representative shape for measurements conducted on this day. On 15 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 ($\sim -\mathrm{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 the size of ice crystals, we do not compare fresh and old snow cases because it is important to note that the date of the old snow case is before fresh snow. This means the studied old snow case is not the aged fresh snow case. Therefore, the change in ice crystal size with its age is not studied in the scope of this paper.

Table 3 Retrieval of the physical characteristics of ice crystals with different shapes in the case of most roughened habits. Numbers in italics indicate minimum root mean square error (RMSE).

Figure 8 Comparison of snow grain size retrieval and best fit of reflectance at three wavelengths of 0.677, 0.873, and 1.649 µm. (a, b, c) Old snow case and (d, e, f) fresh snow case. Effective radius is retrieved at 1.649 µm, and the gray envelope (shading) shows the uncertainty in the retrieved effective radius.

A summary of retrieved effective radii using different ice crystal shapes and corresponding bias and RMSE values is presented in Table 3. The ice crystals with a minimum RMSE value at 1.649 µm are italicized and selected to be used for subsequent calculations of the reflectance factor at 0.677 and 0.873 µm. In Fig. 8, the importance of the 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 it is used to calculate the reflectance factor at 0.677 and 0.873 µm. The results are presented in Fig. 8, with the corresponding RMSE and bias values in the principal plane. The uncertainty of the effective radius retrieval is estimated to be ∼10 %, on the basis of the optimal estimation technique, and is shown by the gray envelope (shading) in Fig. 8. It can be seen that the retrieved effective radius value changes from shape to shape. The difference in retrieved effective radius generally does not exceed 40 %, but in the case of the plate ice crystals, the retrieved effective radius is ∼70 % smaller than the other shapes, e.g., the aggregate of eight columns. This is a significantly large difference. However, these results are presented for the principal plane in which the maximum differences between simulation and measurement are expected. Therefore, the overall bias and RMSE value on all azimuth directions is smaller than what is presented here. It can be seen that the RMSE values at 0.677 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 factor is divided by the measured reflectance.

Figure 9 (a, b, c) The reflectance factor at three wavelengths, namely 0.677, 0.873, and 1.649 µm, from the CAR measurements acquired on 7 April 2008 at Barrow/Utqiaġvik, Alaska, at an altitude of 647 m. (d, e, f) The absolute difference between reflectance of simulation (RSCIATRAN) and that of measurement (RCAR): (RSCIATRAN − RCAR). (g, h, i) The relative difference in percent (%).

7 Comparison of measured and simulated reflectance factor

In this section, we present the results of the comparison of the measured and simulated reflectance factor in the snow–atmosphere system. The simulations, which used the results and findings described in the previous section, were performed by assuming an atmosphere containing O₃, NO₂, and weakly absorbing aerosol, as described in Table 1. The snow layer is assumed to be comprised of an aggregate of eight column ice crystals, with a maximum dimension of 650 µm (effective radius 99 µm) for the case of the old snow, and an aggregate of five plates ice crystals, with a maximum dimension of 725 µm (effective radius 83 µm) for the case of the fresh snow.

Figure 10 The same as Fig. 9 but the measurements by the CAR instrument were performed on 15 April 2008 at 181 m flight altitude over fresh snow.

In Fig. 9, the difference between the simulated and measured reflectance factor 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 $\sim \pm \mathrm{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^∘. Figure 10 is the same plot as Fig. 9 but for fresh snow. The differences between the SCIATRAN simulations and CAR measurements of the reflectance factor are less pronounced in the glint region, compared to those for the old snow.

To assess the accuracy of the simulations over all azimuth angles, the correlation plot and the Pearson correlation coefficient between the measured and modeled reflectance factor are shown in Fig. 11. As it is shown in Fig. 11, the correlation coefficient between reflectance measurements over old snow and simulation is high, ∼0.98. We consider that surface inhomogeneities and related larger shadowing effects at a 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 inhomogeneities are smoothed, and therefore, the old snow case is less affected by surface inhomogeneities). However, in Fig. 11, at 1.649 µm the correlation coefficient is high ∼0.97 for the case of fresh snow, possibly because of there being less sensitivity in this channel to shadowing and atmospheric scattering effects.

Figure 11 The scatterplot with corresponding Pearson correlation coefficient of reflectance factor measured by CAR and simulated by SCIATRAN. (a, b, c) The results for old snow and (d, e, f) fresh snow. Here the color bar represents the number density of the pixels.

8 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 airborne CAR data were acquired by NASA over Elson Lagoon at Barrow/Utqiaġvik, Alaska, during the ARCTAS campaign in spring 2008. The spectral coverage of the airborne measurements is wide (0.3–2.30 µm), comprising important wavelengths relevant for remote sensing applications, such as aerosol retrievals, which could benefit from the results of this study. Measurements obtained at different flight altitudes (∼200, 600, and 1700 m) provide an opportunity to investigate the sensitivity of simulated reflectance to atmospheric parameters.

The SCIATRAN RTM (a phenomenological RTM) was used to simulate the reflectance factor in the snow–atmosphere system and its changes for different snow morphologies (i.e., snow grain size and shape). These simulations explicitly take atmospheric scattering and absorption 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 in order 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 factor, assuming the ice crystals with an aggregate composed of eight column shape agreed well with measurements for the old snow case, has an 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 five plates shape was the most representative ice crystal, with an RMSE value 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, compared to the older snow case at 647 m, are attributed to surface inhomogeneity. The surface inhomogeneity most likely originates from sastrugi. Simulations, in which the snow layer was comprised of ice crystals with a droxtal shape (being semispherical 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 factor, with respect to spectral and directional signatures, produce the measurements well, as evidenced by the high correlation coefficients in the range from 0.88 to 0.98 between measurements (old and fresh snow) and simulations 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 factor from SCIATRAN and CAR measurements is below 0.05. This absolute difference in the off-glint area is smaller in the short-wave infrared 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 satellite or airborne passive remote sensing measurements above snow. To mitigate the relatively large differences between measurements and simulation for the glint condition compared to the 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 the transregional (AC)³ project (Wendisch et al., 2017) that aims to identify and quantify the different parameters involved in the 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 snow–atmosphere system in the algorithms designed to retrieve atmospheric parameters (such as AOT) above polar regions.

Appendix A

For the selected snow models using different ice crystal morphologies, the variation in the snow reflectance R(λ,Ω) at wavelength λ and direction Ω, with respect to the variation δr_e(z) in the effective radius profile r_e(z) along the vertical coordinate z within snow layer, can be presented, neglecting nonlinear terms, by the following equation:

$$\begin{array}{}\text{(A1)}& R\left(\mathit{\lambda },\mathrm{\Omega }\right)=R_{\mathrm{0}}\left(\mathit{\lambda },\mathrm{\Omega }\right)+\underset{\mathrm{0}}{\overset{Zt}{\int }}{W}_{\mathrm{r}}\left(z,\mathit{\lambda },\mathrm{\Omega }\right)\mathit{\delta }{r}_{\mathrm{e}}\left(z\right)\mathrm{d}z,\end{array}$$

where R₀(λ,Ω) and R(λ,Ω) are the reflection functions calculated, assuming an effective radius profile of r_e(z) and r_e(z)+δr_e(z), respectively. The angular variable Ω={θ₀θφ} comprises a set of variables, where θ₀ 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, in the following,

$$\begin{array}{}\text{(A2)}& W_{\mathrm{r}}\left(Z,\mathit{\lambda },\mathrm{\Omega }\right)={\displaystyle \frac{\mathit{\delta }R\left(\mathit{\lambda },\mathrm{\Omega }\right)}{\mathit{\delta }{r}_{\mathrm{e}}\left(z\right)}}\end{array}$$

is the functional derivative of the function R(λ,Ω) with respect to the function r_e(z), which is also called the 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 the 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. (A1) can be used to retrieve the vertical profile of the effective radius within the snow layer in a way similar to that used for the morphology of water droplets (Kokhanovsky and Rozanov, 2012), we restrict ourselves to the assumption of the independence of the altitude r_e. By introducing the weighting function for the absolute variation of the effective radius as follows:

$$\begin{array}{}\text{(A3)}& W_{\mathrm{r}}\left(\mathit{\lambda },\mathrm{\Omega }\right)=\underset{\mathrm{0}}{\overset{Zt}{\int }}{W}_{\mathrm{r}}\left(z,\mathit{\lambda },\mathrm{\Omega }\right)\mathrm{d}z,\end{array}$$

we have the following:

$$\begin{array}{}\text{(A4)}& R\left(\mathit{\lambda },\mathrm{\Omega }\right)=R_{\mathrm{0}}\left(\mathit{\lambda },\mathrm{\Omega }\right)+W_{\mathrm{r}}\left(z,\mathit{\lambda },\mathrm{\Omega }\right)\mathit{\delta }{r}_{\mathrm{e}}.\end{array}$$

The resultant linear relationship is a basic equation used to formulate the 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:

$$\begin{array}{}\text{(A5)}& \mathit{Y}-{\mathit{Y}}_{\mathrm{0}}=\mathbf{K}\left(\mathit{X}-{\mathit{X}}_{\mathrm{0}}\right).\end{array}$$

The components of vectors Y and Y₀ are the measured and simulated reflectance at a 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, and X₀=[ŕ_e] 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:

$$\begin{array}{}\text{(A6)}& \mathrm{\Delta }={∥\mathit{Y}-{\mathit{Y}}_{\mathrm{0}}-\mathbf{K}\left(\mathit{X}-{\mathit{X}}_{\mathrm{0}}\right)∥}^{\mathrm{2}},\end{array}$$

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 as follows:

$$\begin{array}{}\text{(A7)}& \mathit{X}={\mathit{X}}_{\mathrm{0}}+{\left({\mathbf{K}}^T\mathbf{K}\right)}^{-\mathrm{1}}{\mathbf{K}}^T\left(\mathit{Y}-{\mathit{Y}}_{\mathrm{0}}\right).\end{array}$$

In deriving Eq. (7), we have neglected the linearization error which can be significant if X₀ is far from X. To mitigate the impact of the linearization error, we solve the minimization problem given by Eq. (6) iteratively. In particular, instead of Eq. (7), the following is used:

$$\begin{array}{}\text{(A8)}& {\mathit{X}}_n={\mathit{X}}_{n-\mathrm{1}}+{\left({\mathbf{K}}_{n-\mathrm{1}}^T{\mathbf{K}}_{n-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{K}}_{n-\mathrm{1}}^T\left(\mathit{Y}-{\mathit{Y}}_{n-\mathrm{1}}\right),\end{array}$$

where $n=\mathrm{1},\mathrm{2},\mathrm{\dots }$ is the iteration number, K_n−1 and Y_n−1 are the matrix of the 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 the weighting functions and reflectance at flight altitude 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.