Measurements and Modeling of Optical-Equivalent Snow Grain Sizes under Arctic Low-Sun Conditions

The size and shape of snow grains directly impacts the reflection by a snowpack. In this article, different approaches to retrieve the optical-equivalent snow grain size (ropt) or, alternatively, the specific surface area (SSA) using satellite, airborne, and ground-based observations are compared and used to evaluate ICON-ART (ICOsahedral Nonhydrostatic—Aerosols and Reactive Trace gases) simulations. The retrieval methods are based on optical measurements and rely on the ropt-dependent absorption of solar radiation in snow. The measurement data were taken during a three-week campaign that was conducted in the North of Greenland in March/April 2018, such that the retrieval methods and radiation measurements are affected by enhanced uncertainties under these low-Sun conditions. An adjusted airborne retrieval method is applied which uses the albedo at 1700 nm wavelength and combines an atmospheric and snow radiative transfer model to account for the direct-to-global fraction of the solar radiation incident on the snow. From this approach, we achieved a significantly improved uncertainty (<25%) and a reduced effect of atmospheric masking compared to the previous method. Ground-based in situ measurements indicated an increase of ropt of 15 μm within a five-day period after a snowfall event which is small compared to previous observations under similar temperature regimes. ICON-ART captured the observed change of ropt during snowfall events, but systematically overestimated the subsequent snow grain growth by about 100%. Adjusting the growth rate factor to 0.012 μm2 s−1 minimized the difference between model and observations. Satellite-based and airborne retrieval methods showed higher ropt over sea ice (<300 μm) than over land surfaces (<100 μm) which was reduced by data filtering of surface roughness features. ModerateResolution Imaging Spectroradiometer (MODIS) retrievals revealed a large spread within a series of subsequent individual overpasses, indicating their limitations in observing the snow grain size evolution in early spring conditions with low Sun.


Introduction
The enhanced sensitivity of the Arctic climate system regarding global warming, referred to as Arctic Amplification, is associated with several feedback mechanisms [1][2][3]. Numerous with ρ ice representing the density of ice (917 kg m −3 ). For simplification, in the following we use the term snow grain size, which refers to the more accurate term optical-equivalent snow grain size. The snow metamorphism also affects the surface radiative energy budget. More spherical and larger snow grains amplify the absorption of solar radiation and lead to an increase of the surface temperature that in turn accelerates the snow metamorphism. Larger grains allow for a deeper penetration of the incident radiation into the snowpack linked to a higher probability of absorption in the shortwave-infrared (SWIR) spectral range and a decrease of the snow surface albedo, e.g., [11,27].
The SLTSR-based retrieval results were validated against data from seven field-based measurements showing correlation coefficients higher than 0.85 with root mean square errors for r opt and SSA of less than 15 µm and 10 m 2 kg −1 , respectively [39]. A similar correlation coefficient (0.86) was derived from the comparison of the MODIS-based product with ground measurements from six field experiments [40]. Wiebe et al. [40] found maximum differences of 20 µm for undisturbed cases. However, in the presence of cirrus clouds, the retrieved snow grain size tended to be underestimated, while in the case of surface hoar and wind crust an overestimation a maximum difference of 63 µm was observed.
The common basis of many studies is the retrieval of the snow grain size, which is obtained by applying the analytical SSA-snow reflection relationship. This relationship is derived from the asymptotic radiative transfer (ART) approach by Kokhanovsky and Zege [41] assuming a plane surface. The ground-based methods mostly refer to the measurements of the surface albedo, whereas satellite data provide the bidirectional reflectance distribution function (BRDF). However, both reflection properties are influenced by e.g., the surface roughness and the snow grain shape and orientation. Therefore, a comparison of grain size data derived from different observation platforms is significantly affected by the assumption and uncertainty of these parameters.
In general, an increasing surface roughness tends to reduce the surface albedo and leads to a positive bias of the retrieved snow grain size [11,29]. For satellite-based remote sensing of the snow grain size, the deviation of the snow BRDF from that of an ideal plane surface, may lead to an underestimation (overestimation) of the retrieved SSA (r opt ) ranging up to one order of magnitude [42]. The influence of the grain shape on the SSA-albedo/BRDF relationship was explored by several authors, e.g., [27,32,43]. Based on ray tracing simulations, Picard et al. [27] revealed an uncertainty of ±20% of the retrieved SSA from surface albedo measurement when the snow grain shape is unknown. Jin et al. [32] studied the shape effect on satellite-based r opt retrievals and summarized that the directional reflectance is more affected by the grain shape than the albedo, and that the best agreement to measured quantities was found when assuming aggregated snow grains. Assuming a combination of different grain shapes was proposed by Libois et al. [28], since metamorphized snow is mostly composed of a mixture of shapes [27].
The following study compares different methods of ground-based, airborne and satellite-based observations of the snow grain size. The intercomparison is further discussed in relation to modeled data from a numerical weather and climate model, and a parametrization of SSA evolution [44]. In contrast to previous studies on methodical comparisons, e.g., [30], this work applies r opt retrievals on data collected under extreme Arctic conditions (low Sun with SZA about 80°). For these conditions, remote sensing, based on optical measurements, is increasingly challenging because of enhanced measurement and retrieval uncertainties. Since low-Sun observations are prevalent especially in Arctic spring and autumn, and an evaluation of weather and climate models require observations of larger spatial scales, this study estimates the variations of different snow grain size retrievals for these extreme conditions. Section 2 introduces the instrumentation and analyzed data set, which was obtained in the framework of a three-week measurement campaign in the North of Greenland in March/April 2018. Section 3 presents the applied models and retrieval methods to estimate the snow grain size. Section 4 shows the intercomparison with respect to (i) the temporal variability of local snow grain size measurements and modeling, and (ii) the spatial variability based on satellite and airborne observations. Furthermore, the retrieval uncertainties are discussed (Section 5) before a summary and a conclusion are given in Section 6.

PAMARCMiP Campaign
This study is based on measurements performed during the Polar Airborne Measurements and Arctic Regional Climate Model Simulation Project (PAMARCMiP) in 2018. PAMARCMiP 2018 belongs to a series of aircraft campaigns performed within the Arctic region [45] and was conducted together with ground-based observations from 10 March to 8 April 2018. Ground-based and airborne measurements were performed at and in the vicinity of the Villum research station (Station Nord), Greenland (81°36 N, 16°40 W) to document the short-term variability, horizontal and vertical distribution of aerosols and BC in the atmosphere, and concentrations of BC embedded in snow. The airborne activities started on 23 March 2018 and were carried out with the research aircraft Polar 5 [46]. During 14 flights cloud, aerosol [47], and surface properties were quantified by in situ and remote sensing instruments. The observations mainly covered Arctic ocean and the Fram Strait.
Surface properties, as the spectral surface albedo and snow grain size were derived from the spectral modular airborne radiation measurement system (SMART) [48]. The airborne laser scanner RIEGL VQ580 measured the distance to the surface with an accuracy of about 2.5 cm [49]. Out of these data, a 1 × 1 km 2 reference elevation model with a horizontal resolution of 1 m was generated along the flight track. The standard deviation of the relative surface elevation describes the surface roughness. A downward-looking commercial photo camera equipped with a fisheye lens was used to classify the surface conditions. To quantify atmospheric properties, dropsondes of type RD94 [50] were released during the flights. Vaisala HUMICAP humidity and temperature sensors were part of the basis meteorology instrumentation of the Polar 5 aircraft. An airborne Sun photometer with an active tracking system (SPTA) [51] was installed on the top of the aircraft and provided the aerosol optical depth (AOD) at 861 nm and 1026 nm wavelengths. Atmospheric aerosol was also characterized by the Airborne Mobile Aerosol LiDAR (AMALi) system [52] operated in zenith-viewing direction to measure backscatter coefficient profiles at 355 and 532 nm wavelength.

Ground-Based Measurements by the IceCube System
As a ground-based reference, an IceCube instrument was used to derive the SSA of snow over land during PAMARCMiP 2018. The SSA was measured daily at the ground along a fixed 100 m transect located in close vicinity of the Villum research station (distance of 2 km) between 19 March and 4 April, with about 51 samples taken each day. Additionally, broadband surface albedo measurements (300-3600 nm wavelength) were performed by a pair of stationary pyranometers (CM22 by Kipp&Zonen, Delft, The Netherlands) installed close to this IceCube sample line. The manufacturer gives an irradiance uncertainty of about 2%. This uncertainty increases for low-Sun measurements due to the increase of the cosine response error (max. ±3% deviation from ideal at 80°SZA). A second SSA data set was sampled between 22 March and 3 April along a 150 m transect with 5 samples each day about 600 m away from the other transect. With these data sets, temporal and spatial variabilities of the snow grain size and SSA within the course of the campaign were observed.
The IceCube device illuminates a snow sample with a laser diode emitting at 1310 nm wavelength underneath an integrating sphere [38]. A photodiode detects the reflected signal, which is used to calculate the SSA based on radiative transfer simulations with an uncertainty of about 10% for SSA values of up to 60 m 2 kg −1 , which corresponds to a r opt down to 55 µm [38]. The limitations of the IceCube measurement principle for snow samples with smaller grain sizes is related to artefacts, which occur when the snow density is lower than 100 kg m −3 and the radiation may reach the bottom of the snow sample. However, the mean density derived from the weight of the snow sample was 230 ± 30 kg m −3 during PAMARCMiP 2018 and did not fall below this threshold. Recently, Calonne et al. [53] found a systematic factor of 1.3 between SSA derived from IceCube measurements and tomographic images. Optically based SSA retrievals as for the IceCube depend on assumptions about the snow grain shape, such that larger uncertainties than 10% may occur.

Airborne Measurements by SMART
The SMART instrument on board of the research aircraft Polar 5 consists of optical inlets, fibre optics, spectrometers, and a data acquisition system. A set of upward-and downward-looking optical inlets were installed on the aircraft fuselage. The optical inlets were actively stabilized to correct for aircraft movement [48]. The upward and downward spectral radiation was transferred, via optical fibre, from the optical inlets to a set of four spectrometers (two for each hemisphere) covering a spectral range of 0.3 µm to 2.2 µm wavelength with a full width at half maximum of 1-2 and 9-16 nm, respectively [48,54,55]. Radiometric calibrations were performed before and after the field campaign using a NISTcertified (National Institute of Standard and Technology) radiation source (1000 W lamp). In addition, in-field calibrations were applied documenting possible temporal drifts of the SMART sensitivity during the campaign. At large solar zenith angles around 80°to 85°as present during PAMARCMiP, the uncertainty of the measured irradiance at flight level is increased compared to observations performed at smaller SZA. The known components of uncertainty (cosine correction, sensor tilting, absolute calibration, transfer calibration, wavelength accuracy, and dark current subtraction) of SMART were re-evaluated with respect to the large SZAs and the wavelengths applied in the snow grain size retrieval. In particular, the uncertainty of the cosine correction (4%) and the uncertainty of the sensor tilt (2.5%) have a major effect on the overall accuracy of the downward irradiance in the nearinfrared (NIR) wavelength range, since the direct-to-global fraction is approaching unity in this spectral range. Using Gaussian error propagation, the uncertainty of downward and upward irradiance in the NIR summed up to 5.7% and 4.0%, respectively.

Satellite Measurements
Two different approaches for satellite SSA retrievals were considered in the study, based on data from the MODerate Resolution Imaging Spectroradiometers (MODIS) on board of the Terra and Aqua satellites and the Sea and Land Surface Temperature Radiometer (SLSTR) instrument on board of Sentinel-3. SLSTR covers the VIS, NIR, and infrared spectral range with nine spectral channels. The used channels for the snow grain size retrievals (0.55 µm and 1.6 µm) have a spatial resolution of 500 m with a measurement accuracy between 2 and 5% [56]. MODIS obtains data in 36 spectral channels with wavelengths ranging from 0.405 to 14.385 µm. Radiance data (level 1B product MOD02) of three channels (3: 0.47 µm, 2: 0.85 µm, 5: 1.24 µm) are applied for the retrieval. These data have a spatial resolution of 250 m and 500 m, respectively, and show a radiometric accuracy of 1.5% to 3% [40].
For the period and area of the PAMARCMiP observations, SLSTR data were available almost once a day, while from MODIS up to four images per day could be used for the data evaluation.

Measurement Conditions during PAMARCMiP 2018 2.3.1. Sea Ice Conditions
The analysis of aircraft observations focuses on the period 25-27 March 2018, when mostly cloudless conditions prevailed along the flight paths. This restriction to a cloudless period is required to have collocated satellite observations of the surface available. Figure 1 shows the sea ice roughness as derived from the airborne laser scanner along the flight tracks (black lines) for these days.

Meteorological Conditions
The general meteorological situation during PAMARCMiP 2018 was characterized by a high-pressure system over the North Pole and weak lows over North-East Greenland, leading to a period of cloudless conditions in the measurement area between 25 and 27 March 2018. Observations of the sea ice surface properties were performed in different flight altitudes ranging from 50 m to 5 km partly passing the same location. For the radiative transfer simulations in this study the atmospheric conditions were constrained by the airborne observations. Sun photometer data were used to estimate the AOD of the entire atmospheric column. The vertical distribution of the aerosol as indicated by AOD measurements in different altitudes was similar for all three days with a continuous decrease of AOD with flight height. No indication of distinctive aerosol layers up to 5 km was given. To setup the simulations in the NIR spectral range, the AOD was extrapolated to 1100, 1280, and 1700 nm by fitting the Angstrom formula [59] to the measured AOD at 861 and 1026 nm. For determining the columnar AOD, only data from flight sections in the lowest altitude were taken into account. For the three flights, mean columnar AODs at 1700 nm wavelength between 0.01 and 0.03 were obtained, indicating the clean conditions during the three flights. In addition, measurements with the AMALi system did not show any disturbances by clouds or aerosol layers. Only a short sequence on 25 March around 16 UTC was removed from the analysis. The atmospheric profiles of air humidity and temperature were compiled from aircraft and dropsonde data. Dropsondes were released during the flights on 26 and 27 March, while for 25 March the atmospheric profile was derived on basis of the aircraft meteorological sensors during a continuous ascent. The temperature profile on 26 March shows the strongest inversion of all three flights, with −20 • C difference between the surface temperature of −30 • C and the temperature at the inversion height located around 880 hPa, corresponding to an altitude of 1.1 km. The weakest inversion of about 5 K difference was measured on 27 March. The main flight patterns on both days were performed close to 82.5°N over sea ice. On these days, the absolute humidity below 900 hPa pressure level was significantly lower than on 25 March, when the near-surface absolute humidity was affected by areas of open water close to 81.5°N latitude. For this reason, the largest atmospheric impact on the measured radiative quantities due to extinction is expected for the flight on 25 March.

Overview
This section introduces the modeling tools and snow grain size retrieval methods that were applied on satellite and airborne radiation measurements. For better orientation Figure 2 provides an overview of the linkages between the different measurements, retrievals, and models used to estimate the temporal and spatial variability of the snow grain size (or SSA) over sea ice and land surface. Additionally, to the observational results, simulations of the snow grain size metamorphism performed by means of a weather model and a SSA evolution scheme are presented and compared to the temporal evolution of the snow size derived from local IceCube measurements. The retrieval algorithm applied on SMART data is based on atmospheric and snow radiative transfer models (RTM) which were coupled iteratively. The atmospheric RTM provides the direct-to-global fraction ( f dir/glo ) of the solar radiation incident on the snow depending on the atmospheric conditions and the surface albedo. This direct-to-global fraction is set as boundary condition in the snow RTM for creating Look-Up-Tables (LUT) that are used for the snow grain size retrieval.

Snow Radiative Transfer Model-TARTES
The open-source Two-streAm Radiative TransfEr in Snow model (TARTES) [43] was used to simulate the surface albedo in direct and inverse mode. The calculation of the snow albedo from snow grain size in direct mode served as input to atmospheric radiative transfer modeling to assess e.g., the impact of clouds on snow albedo and the impact of the assumed snow grain shape on the retrieved grain size. In inverted mode TARTES was used to retrieve the snow grain size from aircraft albedo measurements.
TARTES simulates the radiative transfer in a snowpack applying the delta-Eddington approximation [60]. The snowpack can be constructed from a predefined number of horizontally homogeneous snow layers defined by their snow density, SSA, and mass fraction of soot. The description of the single-scattering properties of each layer is based on analytical equations given by Kokhanovsky and Zege [41] (see also Section 4).
Libois et al. [43,61] discussed the role of the snow grain shape on the radiative transfer in a snowpack. The grain shape is represented by the absorption enhancement parameter (B), and geometrical asymmetry factor (g G ). B accounts for the photon path length inside the snow grains due to multiple internal reflections, while g G approximates the ratio between forward and backward scattering by the snow grains. Following Libois et al. [43], for particles large compared to the wavelength, the asymmetry factor g can be estimated by: TARTES allows an adjustment of B and g depending on the selected particle shape.
In the simulations of the snow surface albedo, we assumed a single snow layer without soot impurities. The snow grain size was varied between 10 µm and 300 µm. The surface albedo strongly depends on the spectral distribution and the direct-to-global fraction of the incident radiation. This input was provided by an atmospheric radiative transfer model (Section 3.2.2). Generally, for a smooth snow surface the surface albedo increases with increasing SZA due to a higher probability of the photons to be scattered out of the topmost layer of the snowpack at low Sun. Additionally, the forward scattering dominates the asymmetry of scattering, and increases the surface albedo [11]. For low Sun, singlescattering dominates, while for higher Sun the radiation can penetrate deeper into the snowpack corresponding to a higher probability of multiple-scattering. The scattering phase function of the snow particles depends on the snow grain shape. Therefore, the effect of the grain shape on the radiative transfer becomes more relevant for single-scattering than for multi-scattering events, when the angular scattering dependence is increasingly smeared out [11].

Atmospheric Radiative Transfer Model-libRadtran
To calculate the direct-to-global fraction of the incident solar radiation and for the atmospheric correction of the airborne surface albedo measurements, the radiative transfer package libRadtran [62,63] was applied. As a solver for the radiative transfer equation, the Discrete Ordinate Radiative Transfer solver (DISORT) [64] routine was chosen. For the parametrization of the gas absorption, the SBDART model [65] was applied. The extraterrestrial spectrum was taken from Gueymard [66]. Profiles of pressure, temperature, density, and gases were adapted to the airborne observations. The aerosol particle properties were specified by the spectral AOD, derived from Sun photometer measurements, the single-scattering albedo (ω), and the asymmetry factor of the aerosol particles. The latter two parameters were estimated from the Ny-Ålesund AERONET (AErosol RObotic NETwork) data set. We set ω = 0.95 and g = 0.65 as default in the NIR. The impact of the uncertainty of ω and g on the simulated NIR spectral irradiance is low, since the AOD derived for the selected data set did not exceed 0.03 (Section 2). Simulations using an ω of 0.99 and a g of 0.58 resulted in a difference to the default settings of less than 1%.

Weather and Climate Model-ICON-ART
Often, satellite measurements serve as validation of models. However, in terms of snow grain size, there are large uncertainties in both models and remote sensing methods. For this reason, this study compares the results of different observational methods during PAMARCMiP with results from the weather model ICOsahedral Nonhydrostatic model (ICON) [67]. In this way, it is possible to assess whether the model can provide an estimate of snow grain size in the absence of measurements. ICON was developed by the German Weather Service (DWD) and Max Planck Institute for Meteorology (MPI-M). The model system solves the compressible Navier-Stokes equations on an icosahedral grid, which can be seamlessly adjusted in resolution for global and regional simulations. A detailed description of the model can be found in Zängl et al. [67] and Giorgetta et al. [68]. With the extension for Aerosols and Reactive Trace gases (ICON-ART) developed at the Karlsruhe Institute of Technology (KIT), the model can simulate aerosols, trace gases, and related feedbacks [69,70]. The limited area mode, applied here, enables the model to simulate a confined region at high resolution with prescribed lateral boundary conditions. The simulation was run with a horizontal resolution of approximately 3.3 km. The initial state and the boundaries were driven with data from the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF) and fed in at six-hour intervals. ICON currently has two different snow models. The first is a single-layer snow model used for the operational weather forecast. The second is an experimental multi-layer snow model [71], which was applied in this study in a three-layer set up. To investigate the impact of aerosols on the optical properties of snow, the model was extended by the snow grain radius as a new prognostic variable, whereby the aging is based on Essery et al. [72]. In contrast to the original parametrization, ICON-ART uses a lower threshold value discriminating new and aged snow. The applied threshold value is 1 kg m −2 compared to 2.5 kg m −2 .

Parametrization of SSA Evolution
Flanner and Zender [44] parameterized the SSA evolution of dry snow with respect to the effect of the local temperature gradient and the curvature growth following the approach by Legagneux et al. [73]. Based on observational data they proposed an empirical relation of temperature controlled SSA evolution by the fit parameters κ and τ: with t for time and SSA 0 representing the initial SSA at t = 0. Simulations in this study were performed for a set of parameters τ and κ representative for a range of snow temperatures (−37°C to −28°C) and vertical temperature gradients (0 K cm −1 to 0.5 K cm −1 ). The temperature-dependent best-fit parameters for τ and κ were fitted to adapt them to the temperature range during the considered period based on the tabulated data at 0°C, −10°C, −20°C, and −50°C in Flanner and Zenner [44].

XBAER Retrieval of Snow Grain Size Using Satellite-Based Sentinel-3 Data
The eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm is a generic algorithm, which can derive aerosol [74], cloud [75], and surface [76] properties from satellite observations. It has recently been extended to derive snow grain size, snow particle shape, and SSA using the Sea and Land Surface Temperature Radiometer (SLSTR) instrument on board Sentinel-3. The retrieval process is performed using a LUT. In the LUT, snow optical properties are pre-calculated for nine predefined ice crystal particle shapes (aggregate of 8 columns, droxtal, hollow bullet rosette, hollow column, plate, aggregate of 5 plates, aggregate of 10 plates, solid bullet rosette, column) [77]. An atmospheric correction step is applied with a weakly absorbing aerosol type [76] and AOD from Modern-Era Retrospective Analysis for Research and Applications (MERRA) simulation. The aerosol profile is approximated by an exponential function between surface and 3 km altitude. Other trace gas profiles are taken from a monthly latitudedependent mean climatology. Snow grain size and snow particle shape are then obtained by minimizing the differences between theoretical simulations and SLSTR observations of the surface directional reflectances at two wavelengths (0.55 µm and 1.6 µm). The sensitivity study, as presented in Mei et al. [37], shows that the impact of snow particle shape selection on the r opt retrieval is significant, and potential cloud/aerosol contamination introduce an underestimation of r opt . The previous comparison between XBAER derived snow grain size and ground-based measurements of continental snow shows a relative difference of less than 5% [39].

SGSP Retrieval of Snow Grain Size Using Satellite-Based MODIS Data
In this study, the snow grain size and pollution amount (SGSP) retrieval algorithm by Zege et al. [35] was applied to MODIS data. Following Wiebe et al. [40], the SGSP retrieval does not require a priori information on the snow grain shape. Radiances of MODIS (MOD02) measured in three channels (0.47 µm, 0.85 µm, 1.24 µm) are used in this method, which reveals a snow grain size retrieval uncertainty of 10% for SZA lower than 75° [35]. This uncertainty increases up to 20% for SZA = 85° [30].
The SGSP retrieval method uses an analytical asymptotic solution of the radiative transfer equation [41]. Following Zege et al. [35], the black-sky surface albedo α bs (θ 0 ), corresponding to the hemispherical reflectance and assuming only direct illumination, can be calculated as a function of the solar zenith angle θ 0 by: where K 0 represents the escape function determining the angular distribution of radiation, which escapes from a semi-infinite, non-absorbing medium as approximated by Kokhanovsky [78] with: For completely diffuse illumination Equation (4) reduces to: defining the white-sky albedo α ws [35]. According to Kokhanovsky and Zege [41] and Zege et al. [35], y in Equations (4) and (6) can be written as: when considering radiative transfer in a dense snowpack, with χ being the imaginary part of the complex refractive index of ice, wavelength λ, which is taken from Warren and Brandt [79]. A represents the form factor, which depends on the particle shape, and combines the absorption enhancement parameter B and the asymmetry parameter g: Ref. [35] gave a range of A between 5.1 for fractals [80] and 6.5 for spheres. This range of possible values of A contributes to the uncertainty of the retrieved r opt (25%) due to the unknown particle shape. The SGSP retrieval uses an averaged value for A of 5.8 with B = 1.5 and g = 0.84, derived for a mixture of randomly oriented hexagonal plates and columns. To reduce uncertainties using different settings for the satellite retrieval and the TARTES simulations, we set A = 5.8 in both applications.
Since satellites cannot measure the albedo directly to relate the snow albedo and the snow grain size using Equation (4), the SGSP retrieval accounts for the BRDF instead. Satellite-based measurements of the snow surface reflectance are determined by both atmospheric and surface contributions. By considering the atmospheric contribution and assuming the spectral independence of the BRDF, the snow grain size is determined iteratively [35,40]. Further details regarding to the theoretical background of the SGSP retrieval, and the applied equations were given in Zege et al. [35] and Wiebe et al. [40].

Snow Grain Size Retrieval Using Airborne SMART Data
Carlsen et al. [30] applied a modified approach of the SGSP retrieval by Zege et al. [35] to derive the snow grain size from airborne spectral albedo measurements. In their retrieval approach, Carlsen et al. [30] used the spectral albedo ratio (R), which is the ratio between the SMART albedo measurements at λ 1 = 1280 nm and λ 2 = 1100 nm wavelength. Based on Equations (4) and (7) they related the snow grain size to R by: It was argued that using a spectral albedo ratio would reduce the retrieval uncertainty, because wavelength-independent calibration uncertainties of the SMART instrument would cancel each other out [30]. Nevertheless, in this study here also a singlewavelength approach is tested that uses the albedo at 1700 nm wavelength (subsequently called α(1700 nm)-based retrieval).
SMART measures the spectral albedo at flight altitude. As for satellite observations, scattering by atmospheric constituents between surface and aircraft alters the radiation spectrum compared to measurements at surface level. Therefore, an atmospheric correction was applied following the method by Wendisch et al. [81]. It is based on an iterative algorithm, which deployed radiative transfer simulations with the radiative transfer package libRadtran [62,63].
The aircraft measurements by Carlsen et al. [30] were performed over the Antarctic Plateau at high elevation and, thus, in dry air and pristine atmospheric conditions, such that f dir/glo in the NIR spectral range was assumed to be close to unity. This allowed the usage of the black-sky albedo in Equation (9) to retrieve r opt under cloudless conditions [30]. However, for the atmospheric conditions during PAMARCMiP the diffuse incident radiation cannot be neglected, such that the blue-sky albedo (α bs ) needs to be taken into account. The blue-sky albedo can be understood as a linear combination of the black-sky and white-sky albedo: Different to Carlsen et al. [30], this study applies a combination of TARTES and libRadtran simulations to generate LUTs. These LUTs relate blue-sky snow surface albedo and snow grain size for the specific atmospheric conditions (in terms of f dir/glo ) during the PAMARCMiP observations. To estimate r opt , a nonlinear least square method is applied which minimizes the root mean square deviation between the observed and modeled albedo.

Relevance of Atmospheric Effect Correction on SMART Retrieval
Both TARTES and the SGSP retrieval method rely on the same theoretical background based on the formalism by Kokhanovsky and Zege [41]. Figure 3 compares the dependence of snow surface albedo with snow grain size for the different approaches. Neglecting the diffuse incident contribution for the PAMARCMiP conditions, would result in a significant difference of the calculated surface albedo for SZA = 80°and A = 5.8 (Figure 3). For all wavelengths, the parameterized black-sky albedo (dashed lines) using Equation (4) is larger than the results from the TARTES simulations (solid lines) and blue-sky-albedo calculations applying Equation (10) (filled squares), which account for the proper f dir/glo . The directto-global fraction and consequently the offset between the black-sky and blue-sky-albedo are wavelength-dependent, such that R shows also a bias between both methods. This indicates the need for considering the direct-to-global fraction in the retrieval and shows the advantage to use coupled atmospheric and snow radiative transfer models. The atmospheric masking over Arctic snow can contribute to significant uncertainties in the albedo-based r opt retrieval. The atmospheric effects representative for the PAMAR-CMiP conditions are illustrated in Figure 4. The spectral surface albedo was simulated for snow grain sizes between 60 µm and 350 µm (SSA: 9.3 to 55 m 2 kg −1 ) using the TARTES model (gray scaled solid lines in Figure 4). The spectral surface albedo for r opt = 60 µm was set as input for atmospheric radiative transfer simulations with libRadtran to calculate the upward and downward spectral irradiances at 200 m and 3000 m altitude, corresponding to common flight altitudes during PAMARCMiP. The height-dependent albedo calculated from the simulated irradiance spectra are shown as dotted and dashed red lines in Figure 4.
Over bright surfaces, such as snow, the atmospheric masking results in a reduction of the albedo in higher altitudes compared to the surface albedo. In the considered wavelength range, the atmospheric masking is dominated by the extinction of water vapour which is most efficient in the gray marked spectral ranges shown in Figure 4. A minor absorption effect on the albedo spectrum for 3 km flight altitude is still visible outside of these marked areas. Only in the range of the atmospheric window (λ > 1550 nm), gas absorption becomes negligible. For the r opt retrieval wavelengths 1100 nm and 1280 nm (both indicated by a vertical line in Figure 4), the albedo at 3 km altitude shows a reduction of 0.14 and 0.12, respectively, as compared to the default surface albedo. The atmospheric impact on the albedo for 200 m flight altitude is rather small with a bias of −0.01. However, the bias would directly contribute to a r opt retrieval error, if no atmospheric correction was applied. The snow grain size matching with the uncorrected albedo at 1280 nm wavelength at 3 km altitude, for example, would result in an overestimated r opt of about 150 µm (SSA = 22 m 2 kg −1 ) compared to the default 60 µm snow grain size. This clearly demonstrates the relevance of the atmospheric correction when using wavelengths, which are highly affected by water vapour absorption. The uncertainty due to an insufficient atmospheric correction is reduced when applying the α(1700 nm)-based retrieval as shown in Figure 4.

Temporal Variability: Local Observations and Modeling
Daily ground-based snow grain size measurements by the IceCube instrument near the Villum research station were performed during PAMARCMiP over almost three weeks starting on 19 March 2018. At the beginning of the measurement period a hard crust covered with only some millimeter of snow was observed, which resulted from a refreezing period after a massive snow melting event in the end of February 2018. After days of snowfall, a period of dry and mostly cloudless conditions followed, whereby the air temperature did not exceed −25°C. The spatially averaged snow grain size data along line A (100 m, 51 samples) and along line B (150 m, 5 samples) are shown in Figure 5. The error bars indicate the 1-sigma standard deviation calculated from the total set of samples. In general, r opt increased slightly at both sample locations between 44 µm and 72 µm within the three weeks of measurements. The highest variability was observed in the first period of snowfall up to the onset of the cloudless period on 25 March 2018. The day-to-day variation in this first week of observations was stronger than for the following periods, and the spatial variation between the 51 samples along the 100 m transect (line A) covered almost the entire range of r opt -values of the three weeks of measurements. Weak snowfall and blowing snow were reported on 20 March, drifting snow and weak snowfall for 22 March, which might explain the striking variability on these two days. The spatially averaged snow grain size along line B showed mostly higher r opt -values than measured along line A (600 m away), in particular in the first week with snowfall and drifting and blowing snow conditions. For the remaining period, both data sets of line A and line B agreed within the range of the individual standard deviations. The range and the temporal evolution of the measured snow grain size is less strong, with an increase of 15 µm within five days after snowfall, than observed by Carlsen et al. [30] for example. Their measurements on the Antarctic Plateau have shown a more pronounced daily increase in snow grain size after snowfall of about 5.8 µm day −1 (daily SSA decrease: 3.2 m 2 kg −1 day −1 ) under a similar temperature regime (−20°C to −35°C). The snow grain size evolution simulated with ICON-ART is shown in Figure 5 as a solid line. The simulation assumes a growth rate factor of 0.06 µm 2 s −1 for this temperature and snow grain size range as suggested in Essery et al. [72]. The snowfall period before 25 March is well covered by the ICON-ART simulations. However, the growth of the snow particles evolves rapidly in the cloudless period reaching r opt -values up to 110 µm (SSA = 30 m 2 kg −1 ) which is about twice the numbers derived from the IceCube measurements. With the onset of the second short period of snowfall on 30 March, the snow grain size decreased to a value similar to the in situ observations. The comparison shows that the snow grain size of new snow can be well reproduced by ICON-ART. However, the aging process is not well represented by the growth rate factor from Essery et al. [72] for the specific conditions during PAMARCMiP. Therefore, the parametrization of the growth rate factor in ICON-ART was adjusted, such that the simulated snow grain size covers the in situ measurements (red dashed line in Figure 5). For the specific temperature and snow grain size during PAMARCMiP, the original growth parametrization of Essery et al. [72] was applied, but with a reduced growth rate factor of about 0.012 µm 2 s −1 , one fifth of the original value.
In addition to ICON-ART, the parametrization by Flanner and Zender [44] was compared to the observations. For the precipitation-free period starting at the end of 24 March we calculated the snow grain size evolution based on Equation (3) for two scenarios with dT/dz = 0 K cm −1 and dT/dz = 0.5 K cm −1 . From snow pit measurements performed on 24 March, a vertical temperature gradient of about 0.4 K cm −1 was derived. Further snow temperature measurements were conducted at the top of the snowpack covering a range between −28°C and −37°C. The dark gray and blue areas shown in Figure 5 account for this range of snowpack top temperatures, where the upper boundaries of the snow grain size range comprise with the higher snowpack top temperature and the lower boundary with calculations for −37°C. A significant overestimation was observed assuming a vertical temperature gradient of 0.5 K cm −1 (blue shaded areas in Figure 5). The snow grain size for this scenario matches well with the original ICON-ART simulations (black solid line). However, the IceCube measurements show that the vertical gradient effect is less relevant for these low surface temperatures than considered in the parametrization after Flanner and Zender. Using an equilibrium metamorphism (dT/dz = 0 K cm −1 ) would lead to a much better agreement between parameterized and measured snow grain sizes (dark gray area in Figure 5). One of the reasons for this poor representation of snow grain size evolution by the parametrization might be caused by the lower fitting quality of the temperature-dependent parameters, τ and κ from Flanner and Zender [44] to the observed temperature range during PAMARCMiP. For dT/dz = 0 K cm −1 , the original temperature-dependent parameters were described by an exponential decay fitting with a coefficient of determination (R 2 ) larger than 0.99, while for dT/dz = 0.5 K cm −1 the fitting quality was significantly lower [R 2 (τ) = 0.5, R 2 (κ) = 0.99].
The snow metamorphism affects the measured broadband surface albedo, which is shown in Figure 5 Despite the observed bias between modeled and measured snow albedo, we used TARTES simulations to evaluate the change of the broadband surface albedo from period I to period II. The decrease of snow albedo might be caused by the increase of snow grain size and/or the change of the atmospheric conditions. To separate these two effects, the snow albedo for period II was re-calculated assuming a r opt of 50 µm (similar to r opt in period I). As a result, the increase of α by 0.01 indicates a minor effect by the snow grain size variation. For a more detailed investigation of the atmospheric impact on the broadband surface albedo, the snow albedo for period II was re-simulated using TARTES assuming only white-sky albedo (similar to the cloud conditions from period I). The new setup forced the surface broadband albedo to increase by 0.08 which emphasizes the impact of clouds. Consequently, for the discussion on the snow grain size effect on the surface albedo, the atmospheric impact must be separated. The temporal decrease of the surface albedo in Figure 5 was attributed to the cloud impact rather than to the increase of the snow grain size for the PAMARCMiP period and conditions.

Retrieved Maps of Snow Grain Size
Maps of the retrieved snow grain size from the SGSP and XBAER retrieval techniques using MODIS and Sentinel-3 data, as well as the reflectance at 1.24 µm wavelength from MODIS measurements at 11:50 UTC for 25 March are shown in Figure 6. The snow grain sizes estimated from the SMART measurements along the flight track (14-17 UTC) are displayed on each of the panels. They were retrieved from the surface albedo at 1700 nm wavelength. Four MODIS overpasses were evaluated for the period and region of aircraft observations on this day. The different number of valid data points led to an irregular spatial distribution of the snow grain size in each of the four MODIS maps (Figure 6a-d). As illustrated in Figure 6, the main spatial features of the retrieved snow grain size show similar patterns from west to east with lowest r opt -values over land, increasing r opt -values near the eastern coast of Greenland, an area of slightly decreasing r opt (near −9°longitude), and highest values in the most eastern part of the overflown area. Both satellite and airborne observations revealed less variation of the snow grain size over Greenland than over the sea ice. Over Greenland, the retrieved r opt was mostly less than 100 µm (SSA = 33 m 2 kg −1 ) , while r opt over sea ice reached values of up to 300 µm (SSA = 11 m 2 kg −1 ). An exception was found for the map from the MODIS 16:45 UTC overpass, where significant lower r opt -values were retrieved over the sea ice (Figure 6b). At this time, the SZA ranged between 82.4°and 84°for the entire scene. The SZA of the other satellite overpasses were smaller between 79.1°and 81.9°. As discussed earlier, the retrieval uncertainty increases with increasing SZA, which might be one of the reasons for the apparent different spatial snow grain size pattern observed in the late afternoon overpass (Figure 6b). The spatial distribution of the reflectance at 1.24 µm wavelength (Figure 6f), which is sensitive to the snow grain size, shows an increasing surface inhomogeneity in the eastern region with the highest r opt -values. A low reflectance at this wavelength does not necessarily correspond to open water. It might also indicate young ice areas with a possible thin snow layer on top, which causes an overestimation of the derived snow grain size. For example, in the area centred at 81°latitude and −11°longitude, such low reflectances together with high r opt -values were measured, while the AMSR instrument indicates a closed sea ice cover. Furthermore, the measurements might be affected by thin low-level clouds generated above open leads, which were not completely excluded from the data analysis. Limited to the area of the Sentinel-3 overpass, the frequency distributions of r opt are shown for each overpass of MODIS and SLSTR in Figure 7. The 13:50 UTC overpass was excluded in this analysis due to the high amount of unclassified pixel, which would bias the statistics of this case. The plot of the relative frequencies (in r opt -bins of 10 µm) shows two main modes for the three MODIS-based distributions. These two modes represent the lower snow grain sizes over land and the higher numbers retrieved over sea ice. The two morning overpasses revealed similar distributions over sea ice, but some shift of the "land"-mode by 20 µm snow grain size. Corresponding to Figure 6b the relative frequency of the MODIS data from 16:45 UTC revealed the smallest distribution and the smallest r opt -values compared to the other MODIS overpasses. The XBAER retrieval shows a significant smaller variability with a narrower frequency distribution. A narrow mode with a maximum at 120 µm marks the snow grain size derived over land. Over sea ice, there are two further modes (maxima at 140 µm and 180 µm, respectively), with the third mode resulting from the highest r opt -values measured over the most eastern region (Figure 6f). There, the surface is more heterogeneous and indicates an effect of surface roughness on the retrieval results (see Section 5).

Statistical Comparison for Smooth Snow Surfaces
The spatial scales of typical roughness features are below the resolution of the satellitebased observations, which makes it difficult to identify sub-scale roughness features from MODIS or SLSTR data alone. Therefore, observations along the aircraft flight path were used to screen the satellite and SMART data for surface conditions. Areas with increased surface roughness, without snow, or with thin snow layers that are not optically thick, and cases, which were contaminated by low-level clouds were identified by camera observations and laser scanner data. The snow surface class was derived from manually selected red, green, and blue channel thresholds which were set from training samples as applied in Jäkel et al. [82] and Hartmann et al. [83]. Since the laser scanner did not cover the entire flight path, we used also fisheye camera images to estimate the fraction of shadowed and illuminated areas within the individual images as a marker for the roughness of the overflown surface. The ratio of the red and blue channel was calculated for each image pixel. From training images, threshold values of the ratio were defined, which characterized shadowed (ratio < 0.8) and illuminated (ratio > 1.1) pixels. The areal fractions of the shadowed and illuminated pixels ( f sh , f il ) were calculated with respect to the angular resolution of the image pixels [82]. A "smooth surface" was set for cases with f sh < 5% and f il < 5%. The spatial distribution of the remaining SMART data points in Figure 6 shows the filtering result with respect to the surface conditions. In particular, the flight section between −12.5°and 10.0°longitude was identified as a region of rough sea ice.
The statistical measures of the retrieved snow grain size are summarized in Figure 8. For the 25 March, the filtered data were separated into observation above sea ice and land. ICON-Art simulations were also available over the land region. The satellite data were matched to the flight track of the Polar 5 aircraft before the statistical mean, the median, the first and third quartile, and the minimum and maximum values without outliers were calculated. A running average of the SMART measurements was used to account for the spatial resolution of the satellite data. For 25 March (sea ice), the analysis reveals that the interquartile ranges (IQR), indicated by the gray boxes, cover different r opt -ranges, especially for the SGSP retrievals of the MODIS data from 13:50 UTC and 16:45 UTC. The applied SGSP retrieval exhibit no clear bias compared to the other methods as can be concluded from the broad range of retrieved snow grain sizes. The XBAER retrieval shows the smallest IQR, and apart from the 16:45 UTC MODIS overpass, also the lowest mean r opt . For the SMART albedo measurements, both the R-based and the α(1700 nm)-based retrieval method were applied. They revealed differences of the mean r opt -values of 47 µm between both methods for the flight over sea ice, and a 12 µm difference over land. This corresponds to differences in the SSA of about 5 m 2 kg −1 over land and sea ice. Overall, the spread between the mean r opt -values of the different methods is significantly lower over the land surface than over the sea ice ( Figure 8). Apart from the XBAER retrieval, the qualitative differences between the methods are similar for observations over land and sea ice, with lowest (highest) r opt -values for the 16:45 UTC (13:50 UTC) observations by MODIS. The best agreement to IceCube measurements were derived from the 16:45 UTC MODIS overpass where the IQR is within the measurement range of the IceCube data (Figure 8b). All other retrievals and simulations showed a positive bias compared to the in situ measurements. XBAER and the SGSP retrieval of the 13:50 UTC MODIS overpass deviate by a factor of two from the IceCube results. Possible reasons for the deviation of the XBAER results, but also of the SMART retrieval, could originate from the assumption of the snow particle shape when calculating the LUTs. Although for the SMART retrieval, a mixture of grain shapes is assumed, XBAER estimates the snow grain shape in 9 classes [37] in addition to the snow grain size. For the considered area, mostly droxtals were retrieved over the land and the coastal region, and aggregates of 8 columns over sea ice. Implications of the shape effect are further discussed in Section 5. The retrieval results of the MODIS instruments suggest that the snow grain size from the Terra satellite tends to be lower than the r opt -distribution derived from the Aqua satellite. Comparing the MODIS snow grain size within a 2 km radius around the Villum research station with IceCube measurements (45-72 µm, including measurement uncertainty), revealed best agreement with r opt -values from the Terra satellite (53-84 µm) on all three days. In contrast, Aqua showed r opt between 78 to 120 µm for the same period. In all three flights, the α(1700 nm)-based retrieval for SMART revealed smaller values than the R-based method. For 27 March, the ICON-ART simulations covered the entire flight track, because the model-based land mask of ICON-ART classified this near coastal region as land. Similar to the comparison over land (25 March), the model showed low variability (<1 µm standard deviation) and lowest r opt (<100 µm) compared to all other methods. Due to the rather coarse resolution of 3.3 km in this model setup, the small-scale variations present in the observations could not be resolved properly.

Discussion: Implications of Low-Sun Conditions
The different r opt retrieval methods from the satellite and airborne observations are all subject to uncertainties which become increasingly relevant under the low-Sun conditions in the Arctic. After addressing the uncertainties of the SGSP (Section 5.1) and SMART (Section 5.2) retrieval, we are discussing the effect of the choice of snow crystal shape (Section 5.3), and the retrieval wavelength that affect the penetration depth of radiation in the snow (Section 5.4).

Uncertainties of the SGSP Satellite Retrieval
The comparison of the different r opt retrieval results from the satellite-based optical observations show a large spread in particular over sea ice ( Figure 6). The reasons behind this large spread among the SGSP retrieval results of successive overpasses are manifold. First, uncertainties introduced by the measured reflectances can contribute significantly to the retrieval bias. Wiebe et al. have shown that a combination of uncertainties in two or three MODIS channels can result in a snow grain size error of up to 36% for 100 µm grains [40]. In a model study, Zege et al. analyzed the effect of MODIS radiances uncertainties on the SGSP retrieval error. They used a random normally distributed error with the standard deviation of 2% and found a factor of two between reference and retrieved r opt for SZA = 80° [35]. Since the measurement uncertainties are affected by Sun and sensor viewing directions, which were variable between the different MODIS overpasses, we can assume that these uncertainties contribute to the observed variation of the retrieved snow grain size.
The atmospheric correction of both satellited-based retrieval methods does not account for actual atmospheric gas profiles. Instead, mean climatology data are used. As shown in Section 3, the retrieval uncertainty is dependent on the accuracy of the atmospheric correction and increasingly relevant for low-Sun conditions [35]. In contrast to XBAER, which used AOD data from MERRA, the SGSP retrieval assumes climatological-based aerosol profiles. Improper aerosol assumptions may lead to a systematic underestimation of retrieved r opt . They are small (3%) for typical background Arctic aerosol conditions (AOD ≈ 0.05), but increase to up to 37% in the case of Arctic pollution conditions for AOD ≈ 0.11 [37]. Since the AOD did not exceed 0.05 during the three days of observations, the retrieval uncertainty due to aerosol effects is small.
Although the XBAER retrieval obtains r opt and snow particle shape by an iterative optimization approach of measured and pre-calculated reflectances, the SGSP retrieval relies on the radiative transfer theory that relates surface albedo and snow grain size. From satellites, one can only observe the surface from one direction, but a hemispherically integrated value, in terms of the surface albedo, is needed to infer r opt by the SGSP retrieval. This brings forward a problem of inferring surface reflectances in all directions (BRDF) using just one directional observation, and this at three spectral channels used in SGSP. The assumption in the SGSP retrieval is that snow BRDF does not depend on wavelength. This assumption does not hold precisely, as there is a dependence on the real part of the complex refractive index between the used MODIS channels 3, 2 and 5, which proved negligible at smaller SZAs. In the cases where snow BRDF deviates significantly from Lambertian surface (i.e., at higher solar angles), the variation with wavelength becomes significant and propagates into the retrieved grain size. The difference between the higher and lower SZA in accounting for the BRDF in SGSP lies in redistribution of first and higher orders of scattering into the hemisphere with changing wavelength. Therefore, the varying real part of the complex refractive index cannot be considered to be merely a multiplication factor which would be easy to account for.
The larger variability of the retrieved snow grain size distributions among the satellite observations that is shown in Figure 7 is additionally caused by differential macroscopic surface roughness effects on the directional reflectance and surface albedo. Compared to plane surfaces, the directional reflectance is reduced in the forward reflectance peak and enhanced in the backward reflectance [84]. The shading effect in the forward direction is more effective than the reflection to the backward direction that leads to a decrease of the hemispherically integrated albedo [42]. However, the satellite-based retrieval methods rely on a compact and plane snowpack, such that roughness effects mainly observed on the sea ice, lead to an increase of the r opt retrieval uncertainty, which is more pronounced for high SZAs. Larue et al. [18] have shown that the surface reflection is sensitive to the fraction and orientation of the roughness features. It can be assumed that the fraction of these features did not significantly change between the different satellite overpasses. However, the relation between orientation of roughness feature, Sun and sensor viewing direction deviates among the satellite observations, such that shading effects by roughness features differently impact the directional reflectance observed by the satellite sensor.

SMART Measurement Uncertainty and Retrieval Sensitivity
In contrast to the SGSP retrieval, the SMART retrieval algorithm directly applies for the albedo-snow grain size relationship that is derived from the asymptotic radiative transfer theory [41]. Since the albedo represents a hemispherically integrated measure, it is less dependent on alterations of the directional reflection. That includes directional effects by surface roughness on a macroscopic scale, but also the impact of the scattering phase function representing the assumed grain shape on a microscopic scale. Both use of the albedo instead of the directional reflectance for the retrieval and application of the atmospheric correction based on directly measured variables reduce the retrieval uncertainties along with low-Sun conditions compared to satellite-based products.
Two wavelength-dependent approaches were tested which show generally lower r opt values for the α(1700 nm)-based than for the R-based retrieval (Figure 8). The retrieval accuracy is affected by the uncertainty of the surface albedo or the albedo ratio. The total uncertainty of the surface albedo retrieved from the airborne observations is estimated to be about 7.1%. Using the albedo ratio R, the uncertainty reduces to 5.8% as the transition to relative measurements provides independence from the absolute calibration.
To estimate the contribution of the SMART measurement uncertainty on the accuracy of the r opt retrieval we applied combined TARTES and libRadtran simulations to relate snow grain size and surface albedo, exemplary for a SZA of 80°. The simulated surface albedo and the albedo ratio R were biased by the corresponding measurement uncertainties of SMART (±∆α). Finally, the snow grain size was retrieved from the biased surface albedo and R using the predefined LUTs. Figure 9a shows the true (input) snow grain size as well as the retrieved grain size for both directions of the albedo bias (±∆α). The retrieved r opt based on the albedo ratio R reveals a larger deviation from the 1:1 line compared to the retrieval results using α(1700 nm), even though the assumed bias of R (5.8%) is smaller than the bias of α (7.1%). As illustrated in Figure 3, the decrease of the surface albedo at 1700 nm wavelength with increasing r opt is steeper than the decrease of R. Therefore, the effect of the measurement uncertainty is higher for the R-based retrieval. For the studied configuration, this may lead to absolute deviations between true r opt and retrieved r opt which are about three to five times higher than using α(1700 nm) (Figure 9a). The relative deviation (Figure 9b) clearly demonstrates the dependence of snow grain size and measurement uncertainty. For small grain sizes, as those of fresh fallen snow, the retrieved snow grain size could be overestimated by about 100% when applying the Rbased retrieval, while the α(1700 nm)-based retrieval would lead to uncertainties of less than 25% for all considered grain sizes.

Effect of Snow Particle Shape
The snow particle shape directly affects the single-scattering properties in terms of differences in the scattering phase function. The MODIS and SMART-based retrievals of the snow grain size were performed for a mixture of particle shapes using similar settings. In contrast, XBAER retrieves the particle shape simultaneously to the snow grain size. The use of an inappropriate ice crystal shape in XBAER may lead to an error between less 10% to more than 50% in the retrieval of grain size, depending on the particle shape and the grain size value itself [37]. An independent ground-based measurement data set of snow grain shape would certainly be helpful to understand the similarity and diversity between different retrievals. However, these data were not available for PAMARCMiP.
Therefore, the sensitivity of the retrieval methods to the assumed snow particle shape was quantified for the PAMARCMiP specific conditions with a SZA of 80°on basis of TARTES simulations. In TARTES snow albedo was calculated for different shapes such as cylinders, spheroids, cuboids, hexagonal plates with variable aspect ratios (height to the length), and fractals for r opt up to 200 µm. The simulated spectral snow albedo served as input for the snow grain size retrieval which was applied for each shape-specific TARTES simulation. The ratio of the retrieved snow grain size (using the α(1700 nm)-based and the R-based retrieval) and the reference snow grain size (mixed shape) of the TARTES simulations is shown in Figure 10a.
Altogether, for the PAMARCMiP specific conditions, the effect of the unknown snow grain particle shape may lead to uncertainties in the range of ±35% in extreme cases when using LUTs based on calculations for a mixed shape particle type which is higher than reported by Picard et al. [27] with uncertainties of ±20%. The tendency of the deviation strongly correlates with the form factor A. Keeping the surface albedo constant, when absorption is enhanced due to the increase of the form factor A, requires a decrease of the snow grain size to compensate the absorption effect. Furthermore, we can conclude from Figure 10a that there is no clear particle type (e.g., cylinder, spheroid) specific tendency of the snow grain size deviation. Rather, the particle aspect ratio may determine the tendency and magnitude of the snow grain size deviation in the same order than the particle type itself. For example, for hexagonal plates the lowest aspect ratio gave a smaller retrieved snow grain size than the mixed shape approach, while for larger aspect ratios the opposite relation was observed, which clearly is driven by the dependence of the asymmetry parameter g on the particle aspect ratio (see Figure 7 in [26]).
It was found that in general, the relative biases between the α(1700 nm)-based and R-based retrieval methods are almost similar for small and large snow grain sizes. The α(1700 nm)-based retrieval shows only a variability of 4% within the studied range of r opt , while for the R-based retrieval this spread is 2%, as illustrated by the standard deviations of r opt in Figure 10a. Both retrieval approaches show similar results for all shapes, with the R-based retrieval being only slightly higher than the α(1700 nm)-based method. The assumption on the grain shape has a much more critical impact on the retrieval. However, for most retrievals no a priori knowledge of the snow shape is available, and the shape mixture is still the best choice. Therefore, the retrieved r opt should be interpreted also as a shape-equivalent grain size representing the snow albedo that can be calculated assuming a shape mixture. In particular, for satellite-based retrievals, a shape mixture might be the best choice, since the observations cover a large footprint (several hundreds of meters) with natural variability of snow grain shape.
Based on the two extreme snow shapes (hexagonal plates with an aspect ratio of 2 and cylinders with an aspect ratio of 0.25), the retrieval algorithm was adapted to either of both by adjusting the form factor A. These modified retrievals were applied to the case of 25 March 2018 separating observations over land and sea ice. The statistics of the retrieved r opt are given in Figure 10b. Absolute mean differences of about 87 µm (α(1700 nm)-based retrieval) and 115 µm (R-based retrieval) over sea ice were derived, while over Greenland the mean differences decreased to 50 µm and 62 µm, respectively. This promotes the usage of the α(1700 nm)-based retrieval for cloudless conditions, because of its lower sensitivity to the snow grain shape than the R-based retrieval. Figure 10. (a) Ratio of retrieved snow grain size for various particle shapes and the reference mixed shape based on LUTs from TARTES simulations for SZA = 80°. The shape-dependent form factors A are given within the plot next to the r opt -averaged values. The vertical bars indicate the standard deviation of the r opt -averaging. The studied shapes: cylinders (Cyl), spheroids (Sph), hexagonal plates (HexP), cuboids (Cub), and fractals are selected according to the TARTES internal shape list from Libois et al. [43]. The number behind the shape abbreviation gives the aspect ratio of the particle. (b) Retrieved r opt from SMART measurements over sea ice and Greenland on 25 March 2018 assuming the two particle shapes hexagonal plates and cylinders.

Wavelength Choice and Penetration Depth
Using different retrieval wavelengths might result in different r opt estimates, because the penetration depth of the radiation in the snow depends on the wavelength and therefore weights the vertical structure of the snowpack differently. This becomes crucial if the snow layers are stratified, such that a vertical difference in the snow grain size can impose systematic differences in the retrieval. According to the Beer-Lambert law, the radiation decreases exponentially with penetration. The distance in the snowpack where the incident irradiance has decayed to 1/e∼37% of its value is the e-folding depth (z e ). It is used as a measure to quantify, for which layer the retrieved snow grain sizes are representative. For snow, the e-folding depth is calculated by: following Zege et al. [85], with ρ snow and ρ ice representing the densities of ice and snow. The penetration depth increases with decreasing wavelength and snow density, as well as with increasing snow grain size. The two non-absorbing retrieval wavelengths of the SGSP and XBAER algorithm (469 nm and 550 nm) are not sensitive to r opt . They are primarily used to derive the soot concentration (SGSP retrieval) and the snow particle shape (XBAER retrieval), respectively. For the other retrieval wavelength (858-1700 nm) the e-folding depth was calculated for snow densities between 200 kg m −3 and 300 kg m −3 (derived from ground-based snow measurements during PAMARCMiP) and snow grain sizes between 60 µm and 180 µm. For these conditions, the SGSP retrieval refers to snow layers of up to 3 cm depth, while the XBAER and SMART retrieval consider snow layers of less than 1 cm depth. Snow pit measurements of the snow grain size and the snow density in the vicinity of the Villum research station have shown only a low variability (less than 5 µm difference) within the first 10 cm of the snowpack, such that the choice of retrieval wavelength to derive r opt is of minor importance here.

Summary and Conclusions
This study compares snow grain size estimates from different observational methods and models under low-Sun conditions. The analysis is based on airborne and groundbased observations during the PAMARCMiP 2018 campaign hold in the vicinity of the Villum research station, North Greenland, in early spring 2018. The applied methods to retrieve r opt are in general all based on optical measurements making use of the grain size dependent absorption of solar radiation by snow, but in detail depend on the specific instrument, which cover ground-based in situ (IceCube), airborne (SMART) and satellite observation (MODIS on Aqua and Terra, SLSTR on Sentinel-3). The different retrieval methods rely on the asymptotic radiative transfer theory [41] applied on airborne albedo and MODIS reflectance measurements (SGSP retrieval) [35], as well as a minimizing approach of measured SLSTR and pre-calculated reflectances for variable grain sizes and shapes (XBAER retrieval).
The snow grain size retrieval of the airborne SMART instrument accounts for the directto-global fraction of the solar radiation incident on the snow by coupling an atmospheric and a snow radiative transfer model. The retrieval was applied for two wavelength settings, (i) an albedo ratio method, and (ii) a new single-wavelength approach using the albedo at 1700 nm wavelength. The reduction of the retrieval uncertainty promotes the usage of this single-wavelength retrieval approach in combination with the coupled atmosphere and snow model.
Moreover, the locally measured r opt evolution was compared to r opt simulations from the ICON-ART model and a parametrization proposed by Flanner and Zender [44]. To our knowledge, these different methods have not been compared at high latitudes (low-Sun conditions) before. In particular, the retrievals using albedo and reflectance measurements are subject to significant uncertainties due to the large SZA of about 80°a s present during the PAMARACMiP campaign. However, conditions with low Sun are common in early spring in the central Arctic. Therefore, this comparison of different approaches demonstrates the consequences of retrieval uncertainties for evaluating the snow evolution. Local in situ measurements over the three-week period of the PAMARCMiP campaign revealed a minor increase of r opt compared to previous measurements on the Antarctic Plateau [30] under a similar temperature regime. The r opt evolution modeled by ICON-ART showed good performance for the time frame of snowfall events. In the cloudless period of the campaign, in contrast to the IceCube in situ data, the modeled r opt doubled its value within five days. Adjusting the growth rate factor to 0.012 µm 2 s −1 led to the best agreement with the in situ data. Additionally, the parametrization after Flanner and Zender [44] showed a significant overestimation of the r opt evolution when assuming a vertical temperature gradient close to the measured gradient of about 0.4 K cm −1 . This indicates certain weaknesses caused by the limited derivation of the best-fit-parameters κ and τ, or the poor representation of the curvature growth for these low temperatures (T < −28°C).
Three days of cloudless conditions were selected to compare ground-based, aircraft and satellite observations of r opt . Measurement flights over the Fram Strait performed on 25 March 2018, indicated higher and more variable r opt -values over the sea ice (r opt < 300 µm) than over land (r opt < 100 µm), which was also deduced from the two satellite-based retrievals, XBAER (SLSTR on Sentinel) and SGSP (MODIS on Aqua and Terra). The statistical analysis of the filtered satellite data covering the flight path of the Polar 5 aircraft over smooth snow surfaces showed mean r opt differences up to 100 µm between the successive overpasses. For land surface measurements near the Villum research station, snow grain size from the Terra satellite (r opt : 53-84 µm) showed a better agreement to the ground-based IceCube data set (r opt : 45-72 µm) than the Aqua product (r opt : 78-120 µm). The difference between XBAER and SGSP snow grain size is larger compared to the difference between SGSP and SMART retrieval, probably due to the assumption of the ice crystal shape. Both SMART retrieval approaches deviated by up to 40% from each other, but ranged between the MODIS derived extremes with better agreement of the α(1700 nm)-based retrieval with the IceCube measurements.
Filtering of the data with respect to smooth sea ice surface conditions did not necessarily improve the comparison between the different retrievals and observations. Measurement uncertainties at low-Sun conditions and the fact that successive satellite overpasses are taken under different Sun and observation geometries, make an additional contribution to the large spread of satellite results. This shows their limitations in studying the day-today evolution of the snow grain size under low-Sun conditions in particular over sea ice. As shown here for one case of PAMARCMiP, the differences of retrieved r opt between two overpasses exceeds the typical evolution of snow grain size by snow metamorphism.
Potential retrieval uncertainties based on the airborne SMART observations were analyzed. The findings of this analysis may serve as recommendations also for satellitebased applications. We propose (i) to apply an atmospheric correction, (ii) to calculate LUTs of the blue-sky albedo, instead of assuming a black-sky albedo, (iii) to consider roughness features and their spatial proportion by collocated laser scanner and/or imaging methods covering a similar FOV, (iv) to make use of suitable wavelengths in the SWIR to use the strongest sensitivity on r opt and lower dependence on atmospheric extinction, and (v) to use a form factor representing a mixed-type of grain shapes.