The retrieval of snow properties from SLSTR Sentinel-3 – Part 1: Method description and sensitivity study

The eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm has been designed for the top-of-atmosphere reflectance measured by the Sea and Land Surface Temperature Radiometer (SLSTR) instrument on board Sentinel-3 to derive snow properties: snow grain size (SGS), snow particle shape (SPS) and specific surface area (SSA) under cloud-free conditions. This is the first part of the paper, to describe the retrieval method and the sensitivity study. Nine pre-defined SPSs (aggregate of 8 columns, droxtal, hollow bullet rosette, hollow column, plate, aggregate of 5 plates, aggregate of 10 plates, solid bullet rosette, column) are used to describe the snow optical properties. The optimal SGS and SPS are estimated iteratively utilizing a look-up-table (LUT) approach. The SSA is then calculated using another pre-calculated LUT for the retrieved SGS and SPS. The optical properties (e.g., phase function) of the ice crystals can reproduce the wavelength-dependent and angular-dependent snow reflectance features, compared to laboratory measurements. A comprehensive study to understand the impact of aerosols, SPS, ice crystal surface roughness, cloud contamination, instrument spectral response function, the snow habit mixture model and snow vertical inhomogeneity in the retrieval accuracy of snow properties has been performed based on SCIATRAN radiative transfer simulations. The main findings are (1) snow angular and spectral reflectance features can be described by the predefined ice crystal properties only when both SGS and SPS can be optimally and iteratively obtained; (2) the impact of ice crystal surface roughness on the retrieval results is minor; (3) SGS and SSA show an inverse linear relationship; (4) the retrieval of SSA assuming a non-convex particle shape, compared to a convex particle shape (e.g., sphere), shows larger retrieval results; (5) aerosol/cloud contamination due to unperfected atmospheric correction and cloud screening introduces underestimation of SGS, “inaccurate” SPS and overestimation of SSA; (6) the impact of the instrument spectral response function introduces an overestimation into retrieved SGS, introduces an underestimation into retrieved SSA and has no impact on retrieved SPS; and (7) the investigation, by taking an ice crystal particle size distribution and habit mixture into account, reveals that XBAER-retrieved SGS agrees better with the mean size, rather than with the mode size, for a given particle size distribution.

1 Introduction

Snow properties such as snow albedo, snow grain size (SGS), snow particle shape (SPS), specific surface area (SSA) and snow purity (Warren and Wiscombe, 1980; Painter et al., 2003; Hansen and Nazarenko, 2004; Taillandier et al., 2007; Gallet et al., 2009; Battaglia et al., 2010; Gardner and Sharp, 2010; Domine et al., 2011; Liu et al., 2012; Qu et al., 2015; Baker, 2019; Pohl et al., 2020a) show large variabilities temporally and spatially (Kukla et al., 1986). They play important roles in the global radiation budget, which is critical to some well-known phenomena such as Arctic amplification (Serreze and Francis, 2006; Domine et al., 2019). Satellites offer an effective way to understand the surface–atmosphere processes and corresponding feedback mechanisms on the regional, continental and/or global scales (Konig et al., 2001; Pope et al., 2014). Satellite-derived snow products (e.g., SGS, SPS and SSA) are particularly important for short-term hydrological, meteorological and climatological modeling (Livneh et al., 2009). A high-quality snow property data product can also be applied to derive aerosol optical thickness (AOT) over the cryosphere (Mei et al., 2020a). High-quality satellite-derived snow products and their by-products are also important for the creation of long-term climate data records (SSMC, 2014), which enable better investigation and interpretation concerning global climate change (Konig et al., 2001). However, both the definition and the corresponding data accuracy of SGS are poor (Langlois et al., 2020), and there is no existing SPS satellite product. The lack of good information on SGS and SPS leads to a low quality of SSA (Gallet et al., 2009). The accuracy of SGS, SPS and SSA limits the model performance for the prediction of snow properties related to climate change issues. Lack of information on SGS and SPS also restricts the accuracy of snow bidirectional reflectance estimation, which further limits the retrieval possibilities of aerosol and cloud properties above snow (Mei et al., 2020a, b).

A comprehensive overview of remote sensing of SGS, SPS and SSA can be found in many previous publications (e.g., Li et al., 2001; Stamnes et al., 2007; Koren, 2009; Lyapustin et al., 2009; Dietz et al., 2012; Wiebe et al., 2013; Frei et al., 2012; Mary et al., 2013; Kokhanovsky, et al., 2019; Xiong and Shi, 2018). The variation in SGS leads to the large variability in top-of-atmosphere (TOA) reflectance in near-infrared (NIR)/shortwave infrared (SWIR) spectral ranges, and SPS shows a strong impact on TOA reflectance in visible channels (Warren and Wiscombe, 1980). Different retrieval algorithms have been developed for different instruments. For instance, the MODIS Snow-Covered Area and Grain size (MODSCAG) retrieval algorithm and Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm have been used to derive SGS using MODIS and VIIRS instruments (Painter et al., 2003, 2009; Lyapustin et al., 2009).

Snow particle shape is another important parameter which affects the estimation of snow properties, such as albedo (Räisänen et al., 2017; Flanner and Zender, 2006), because ice crystals with different shapes have different optical properties (Jin et al., 2008; Yang et al., 2013). The absorption and extinction cross sections of an ice crystal can be described as a function of size, shape, and the refractive index at a given wavelength (van de Hulst, 1981; Mischenko et al., 2002, and references therein). Natural snow consists of grains, depending on temperature, humidity and meteorological conditions, which have numerous different shapes (Nakaya, 1954). SPSs have been classified into different categories; the classification has been increased from 21 (Nakaya and Sekido, 1938) to 121 (Kikuchi et al., 2013) categories. Although a spherical-shape assumption is typically used for field measurements (Flanner and Zender, 2006; Donahue et al., 2020), this approximation is not recommended to be used in retrieval algorithms of satellite measurements because it leads to large differences between observed and simulated wavelength-dependent snow bidirectional reflectance, especially at visible wavelengths (Leroux and Fily, 1998; Aoki et al., 2000; Jin et al., 2008; Dumont et al., 2010; Libois et al., 2013). Improper wavelength-dependent snow bidirectional reflectance caused by a predefined SPS leads to low-quality satellite retrieval results. Some attempts to derive SPS in the ice cloud can be found in previous publications (McFarlane et al., 2005; Cole et al., 2014).

According to Legagneux et al. (2002), SSA is defined as the surface area of ice crystal per unit mass; i.e., $\text{SSA}=A_{\mathrm{t}}/\mathit{\rho }V$, where A_t and V are total surface area and volume, respectively, and ρ is the ice density. SSA includes information on both SGS and SPS, and it is often used to describe the surface area available for chemical processes (Taillandier et al., 2007; Domine et al., 2011; Yamaguchi et al., 2019). SSA is reported to have a good relationship with snow spectral albedo at the shortwave infrared wavelengths (Domine et al., 2007). Optical methods are routinely used to measure SSA in the field (Gallet et al., 2009). Empirical equations have been proposed to describe the change in SSA (Legagneux and Domine, 2005; Taillandier et al., 2007). Few attempts have been made to derive SSA from satellite observations (Mary et al., 2013; Xiong et al., 2018).

This paper presents a new retrieval algorithm to derive SGS, SPS and SSA from satellite observations. In a snow–atmosphere system, satellite-observed TOA reflectances are affected by numerous snow and atmospheric parameters. The parameters, which will be estimated in the framework of the eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm, will be called the target parameters. Other parameters, which the TOA reflectance also depends on, will be called the model parameters. In the case of the XBAER algorithm, the target parameters are SGS, SPS and SSA, whereas the model parameters are aerosol loading, cloud optical thickness and gaseous absorption. Throughout the paper, SGS will be characterized by an effective radius. Following Baum et al. (2011), the effective radius is defined as $\mathrm{3}V/\left(\mathrm{4}{A}_{\mathrm{p}}\right)$, where V and A_p are the volume and average projected area, respectively. As can be seen in the case of a spherical particle, the effective radius is equal to the radius of the sphere. The general concept of the retrieval algorithm is to use simultaneously spectral and angular reflectance measurements, which are sensitive to SGS and SPS. The spectral channels used in the XBAER algorithm are 0.55 and 1.6 µm. Both nadir and oblique observation directions from the Sea and Land Surface Temperature Radiometer (SLSTR) are used. An optimal SGS and SPS pair is achieved by minimizing the difference between measured and simulated atmosphere-corrected surface reflectances. SSA is then calculated based on the retrieved SGS and SPS. Nine predefined SPSs (aggregate of 8 columns, droxtal, hollow bullet rosette, hollow column, plate, aggregate of 5 plates, aggregate of 10 plates, solid bullet rosette, column) (Yang et al., 2013; see Table 1) are used to describe the snow optical properties and to simulate the snow surface reflectance at 0.55 and 1.6 µm at two observation angles.

Table 1 Snow particle shape provided in Yang et al. (2013) database. The abbreviations introduced here will be used later.

There are three points we would like to emphasize to avoid misunderstandings between the snow science community and remote sensing community.

• Usage of the Yang et al. (2013) database for ice crystals in the air (ice cloud) and on the ground (snow). The optical properties of ice crystals presented by Yang et al. (2013) have been widely used to study ice clouds. In recent publications, it has been demonstrated that they can also be used for snow studies (Räisänen et al., 2015; Pirazzini et al., 2015; Saito et al., 2019; Schneider et al., 2019; Pohl et al., 2020b). In fact, the single-scattering properties of ice crystals in the Yang et al. (2013) database are determined solely by the given particle size, shape and refractive index. They can be used to describe the optical properties of both snow particles and ice cloud particles when the particle models represent the aforementioned optical/physical properties (Saito et al., 2019; Masanori Saito, personal communication, 2021).

• Snow particle shape observed from field measurements and derived from satellite observations. For scientists working in a laboratory or on campaign-based studies, the best way to obtain an image of snow is to use an X-ray microtomography or confocal scanning optical microscope/scanning electron microscope (Hagenmuller et al., 2016; Baker et al., 2019; Ian Baker, personal communication, 2021). In a field measurement and its related application areas (e.g., calculation of snow albedo), a spherical-shape assumption is widely used because it is easier to derive other snow properties such as SSAs and snow albedo based on this assumption, compared to using other more complicated shapes (see Appendix). The assumptions of a spherical and non-spherical shape have much less impact on the estimation of snow albedo compared to the bidirectional reflection features of snow (Grenfel and Warren, 1999; Dumont et al., 2010), because SPS has a significant impact on the ice crystal phase function while it has a relatively weak impact on the snow extinction/absorption coefficient (Jin et al., 2008). However, the spherical shape cannot be used to provide typical bidirectional reflection features of snow with the required accuracy (Jin et al., 2008; Dumon et al., 2010; Jiao et al., 2019), which is the fundamental basis for deriving snow properties from satellite remote sensing techniques. Thus, more complicated SPSs, such as those proposed by Yang et al. (2013), are recommended to use in the simulations of the angular distribution of snow reflectance. Besides, both snow albedo and directional reflectance are affected by other factors such as how single particles aggregate.

• SGS and SSA. Although what the definition of a snow grain constitutes is an ongoing debate in different communities, SGS and SPS are two fundamental inputs for any radiative transfer model, which is the basis for the satellite retrievals (Langlois et al., 2020). Typically, the SSA is preferable within the snow science community because SSA is commonly used in further applications based on field measurements. We note, however, according to the definition of SSA, for a given SPS, a unique relationship between SGS and SSA can be derived. SPS is the intermediate but fundamental parameter needed to retrieve SSA in our XBAER algorithm.

This paper is structured as follows: observation characteristics of SLSTR and the laboratory measurements used for sensitivity studies are described in Sect. 2. The theoretical background and the ice crystal database (Yang et al., 2013) are presented in Sect. 3. Section 4 describes the eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) algorithm. The results of a comprehensive sensitivity study using SCIATRAN (Rozanov et al., 2014) simulations are presented in Sect. 5. The conclusions are given in Sect. 6.

2 Data

2.1 SLSTR instrument

The satellite data will be used in two ways throughout the two companion papers (this paper and Mei et al., 2021). In the first part, we perform a statistical analysis of the SLSTR observation/illumination geometries to select realistic settings for the sensitivity study. In the companion paper, the satellite measurements will be used as the inputs of the XBAER algorithm to derive the research satellite products of SGS, SPS and SSA.

The SLSTR instrument on board the European Space Agency (ESA) satellite Sentinel-3 is the successor of the Advanced Along-Track Scanning Radiometer (AATSR) instrument and is used to maintain continuity with the (A)ATSR series of instruments. SLSTR has the heritage of AATSR instrument characteristics, especially the dual-viewing observation capabilities and wavelength settings. In order to have a reasonable setting for observation/illumination geometries in the sensitivity study, we perform a statistical analysis of the SLSTR observation geometries (solar zenith angle, SZA; viewing zenith angle, VZA; relative azimuth angle, RAA), similarly to Mei et al. (2020a). This analysis is essential because (1) it provides a realistic setting of observation/illumination geometries in our sensitivity studies and (2) it helps us to have a complete understanding of the observation-/illumination-related surface/atmospheric properties. Here the definition of RAA has been harmonized with SCIATRAN (Rozanov et al., 2014); namely, the RAA value is equal to 0^∘ under a strict glint condition. The statistical analysis has been performed using observations made over Greenland during April and September 2017. April and September are reported to be representative months of the Arctic (Mei et al., 2020a). Please note that these two months are picked to represent the SLSTR observation characteristic with a typical solar illumination angle; the change in underlying surface properties plays no role in such a selection. Figure 1 shows the frequency of SLSTR observation geometries. The upper panels show the SZA with SLSTR nadir and oblique observations for April and September. We can see that the SZA occurs frequently with a value of 70^∘ for selected months. The VZA and RAA for oblique observation mode are typically around 55^∘ and in a range of [110^∘, 170^∘], respectively. The observation geometries for nadir observation show relatively large variabilities due to larger swath width compared to for oblique observations (1400 vs. 700 km). Larger SZAs can be found especially at the edge of the swath. The VZA and RAA for oblique observation mode are typically in ranges of [0^∘, 55^∘] and [70^∘, 140^∘], respectively. According to the statistical analysis, a combination of SZA, VZA and RAA of 70, 30 and 135^∘ for nadir observation and 70, 55 and 135^∘ for oblique observation can be a reasonable setting for the SLSTR observation geometries for the sensitivity study.

Figure 1 Upper panels are the histograms of SZA for SLSTR observations: (a) nadir during April, (b) oblique during April, (c) nadir during September, and (d) oblique during September. Lower panels are the polar plots of (VZA, RAA) probability for AATSR observations: (e) nadir during April, (f) oblique during April, (g) nadir during September, and (h) oblique during September.

2.2 Laboratory measurements

Laboratory measurements of the bidirectional reflectance of snow samples contain important information about the dependence of the angular structure of snow reflection on the lighting geometry, wavelength and snow physical properties. The comparison of measured and modeled bidirectional reflectance helps to establish the conceptual ideas for the retrieval algorithm. For this comparison, we have selected measurements of fresh and aged snow samples presented by Dumont et al. (2010) and Peltoniemi et al. (2009), respectively.

The fresh snow sample, a cylinder of 30 cm diameter and 12 cm height, was taken from a new wet snow layer at Col de Porte (Chartreuse, France) at 1300 m above sea level during January 2008 (Dumont et al., 2010). The sample was stored in a cold room at −10 ^∘C for 1 week to avoid metamorphic effects during the ensuing measurements. To obtain the bidirectional reflectance factor (BRF), the snow sample was illuminated by a monochromatic light source at an incidence zenith angle of 60^∘. The spectral BRF between 500 and 2600 nm was measured at viewing zenith angles of 0, 30, 60, and 70^∘ and relative azimuth angles 0, 45, 90, 135, and 180^∘ by a spectrogonio radiometer developed at the Laboratoire de Planétologie de Grenoble, France, and using a Spectralon® and an Infragold® sample as a reference (see Dumont et al., 2010, for further details).

The aged snow sample, a cuboid of more than 10 cm height, was taken from an old dry snow layer at Masala, Finland, and brought into a warm laboratory. The spectral BRF between 350 and 2500 nm was measured during the aged process by the Finnish Geodetic Institute field goniospectropolariphotometer (FIGIFIGO) and using a Labsphere Spectralon® 99 % white reference plate. For illumination, a 1000 W Oriel Research quartz tungsten halogen lamp at a zenith angle of 60^∘ was utilized (Peltoniemi et al., 2009). The spectral BRF was obtained at viewing zenith angles of up to 70^∘ in 1^∘ resolution and at relative azimuth angles of 0, 90, 130, 160, 180, 270, 310 and 340^∘. The first and last measurements were made in the principal plane, indicating minor metamorphism in the snow layer during the measurement.

3 Dependence of snow reflectance on target parameters

A comprehensive data library (Yang et al., 2013) containing the scattering, absorption and polarization properties of ice particles in the spectral range from 0.2 to 15 µm was used to calculate radiative transfer through a snow layer (Pohl et al., 2020b). A full set of single-scattering properties is available for nine ice crystal habits presented in Table 1. The maximum dimension of each habit ranges from 2 to 10 000 µm in 189 discrete sizes.

The optical properties of ice crystals depend on wavelength, ice crystal size and shape. Maximal dependence of the single-scattering albedo on the particle size is observed in the spectral ranges where ice absorption cannot be neglected. The asymmetry factor depends on the particle size for the whole spectral range. This dependence can be weaker or stronger at a selected wavelength depending on SPS (see Yang et al., 2013, for details).

To better illustrate the impact of SGS and SPS on the radiative transfer through a snow layer, we have calculated the reflectance of the snow layer consisting of droxtals, aggregates of 8 columns, hollow columns and plates with crystal surface roughness conditions as severely roughened. The simulations of snow reflectance were performed using the radiative transfer package SCIATRAN (Rozanov et al., 2014). The snow layer was defined as a layer directly over a black surface, with snow optical thickness of 500 and a snow geometrical thickness of 1 m. The snow layer is assumed to be vertically and horizontally homogeneous without any surface roughness and composed of monodisperse ice crystals. The impact of snow impurities and scattering/absorption processes in the atmosphere was neglected at this stage. The reflectance of the snow layer as a function of the effective radius of ice crystals at wavelengths of 0.55 and 1.6 µm is presented in Fig. 2. The calculations were performed for typical SLSTR instrument observation/illumination geometries (see Sect. 2.1), with SZA, VZA and RAA equal to 70, 30 and 135^∘ (scattering angle 129^∘).

Figure 2 Reflectance of snow layer at 0.55 and 1.6 µm calculated assuming different SPSs. Observation/illumination geometry: SZA, VZA and RAA were set to 70, 30 and 135^∘, respectively.

There are a couple of criteria we considered for the selection of the optimal wavelengths (0.55 and 1.6 µm) in the XBAER algorithm, for the purpose of creating a long-term satellite snow property dataset with good and stable accuracy.

• We took the overlap channels between AATSR and SLSTR because a consistent long-term satellite snow dataset is possible only when the same algorithm can be applied to both AATSR and SLSTR instruments. In particular, the overlap channels between AATSR and SLSTR are 0.55, 0.66, 0.87, 1.6, 3.7, 10.85 and 12 µm.

• Picking up wavelengths for which the contribution of thermal emission can be ignored, then 0.55, 0.66, 0.87 and 1.6 µm remain.

• Deleting the channel 0.66 µm to avoid the potential impact of O₃ absorption, after that, 0.55, 0.87 and 1.6 µm remain.

• We take into account that the retrieval algorithm is a two-stage algorithm; namely, first it uses channels with minimum impact of the ice crystal shape to retrieve the grain size, and then it selects the shape using channels with minimum impact of grain size. Accounting for the fact that the 0.87 µm channel is impacted by both size and shape, 0.55 and 1.6 µm channels were picked up for the retrieval.

The right panel of Fig. 2 demonstrates the strong dependence of the snow layer reflectance at 1.6 µm on the SGS. One can also see that the dependence of snow reflectance on SPS cannot be neglected. In particular, the same reflectance can be obtained with a combination of different SGS and SPS. For instance, one can see from the right panel of Fig. 2 that the reflectance of the snow layer consisting of droxtals with SGS=200 µm or of plates with SGS=65 µm equals ∼0.035 in both cases. Thus, assuming different SPSs, the values of retrieved SGS can differ three times. The left panel of Fig. 2 demonstrates the dependence of the snow layer reflectance at 0.55 µm on SGS and SPS. It can be seen that the dependence of reflectance on SGS is very weak for droxtals and aggregates of 8 columns. However, reflectance at 0.55 µm decreases with an increase in SGS for hollow columns and plates. The weak oscillations for the reflectances at 0.55 µm can be explained by the joint impact of oscillations in the single-scattering albedo and elements of the scattering matrix presented in the original database. Although the reason for the oscillation in the database is unclear, it is unlikely due to physical phenomena (Masanori Saito, personal communication, 2021).

To illustrate this point, the dependence of the phase function at a 129^∘ scattering angle on SGS is shown in the left panel of Fig. 3. The phase functions (F11 element of the scattering matrix) were extracted from the original database. According to the left panel of Fig. 3, the dependence of snow surface reflectance at 0.55 µm on SGS and SPS is caused mainly by the phase function of ice crystals. Weak oscillations can also be found.

Figure 3 Left panel: phase function at 0.55 µm for scattering angle of 129^∘, extracted from the original database (Yang et al., 2013) as a function of effective radius. Right panel: reflectance of snow layer at 1.6 µm consisting of aggregates of 8 columns, calculated assuming that (1) phase function depends on the effective radius (black line), (2) phase function is constant corresponding to the effective radius of 15 µm (red line), (3) same as (2) but for effective radius of 150 µm (green line) and (4) same as (2) but for effective radius of 1150 µm (blue line).

The above analysis shows that accurate retrieval of SGS requires adequate information about SPS and accounting for the dependence of the phase function on SGS. To better illustrate the impacts of SGS on the ice crystal phase function, we calculated reflectance at 1.6 µm with different SGS values. The right panel of Fig. 3 represents the reflectance of the snow layer, consisting of aggregates of 8 columns, calculated accounting for the dependence of the phase function on the effective radius (black line) and assuming a constant phase function for three selected effective radii equal to 15, 150 and 1150 µm (red, green and blue lines, respectively). It can be seen that the accurate simulation of snow reflection requires accounting for the dependence of the phase function on SGS.

The main findings of the presented investigations can be formulated as follows:

• Reflectance of a snow layer depends on both SGS and SPS.

• Accurate simulation of snow surface reflectance requires accounting for the dependence of the phase function on SGS.

• Spectral channels in the visible spectral range are more sensitive to SPS compared to SGS.

• Spectral channels in the near-infrared spectral range are more sensitive to SGS compared to SPS.

Although the global classification of snow crystal, ice crystal, and solid-precipitation particles suggested in Kikuchi et al. (2013) consists of 121 particle types, we restrict ourselves, in the retrieval algorithm, to nine shapes of ice crystals, for which optical characteristics are represented in the database (Yang et al., 2013). And these nine shapes have been proven, especially from satellite observations, to be able to reproduce typical wavelength/angular features of snow reflectance in reality (Räisänen et al., 2015; Pirazzini et al., 2015; Saito et al., 2019; Schneider et al., 2019; Pohl et al., 2020b). To further illustrate that the selected dataset is able to reproduce the BRF of different snow types, we compared the simulated and measured BRF of fresh (Dumont et al., 2010) and aged (Peltoniemi et al., 2009) snow samples. To reproduce the spectral BRF by SCIATRAN, we use the setup described above in this section and adjust the SGS for each SPS by minimizing the deviation between simulated and measured reflectance at 1.6 µm. Figure 4 shows the simulated BRF in the principal plane at 0.55 µm of fresh and aged snow samples, as well as the respective measurements. The BRF is defined as $\mathit{\pi }I/F$, where I is the reflected radiance and F is the incident irradiance. According to Fig. 4a, for fresh snow, plates are the best shape to reproduce the measured BRF in the vicinity of the forward scattering peak but plates underestimate the BRF at higher viewing zenith angles in the backscattering region. Here, shapes of hollow bullet rosette, hollow column and aggregate of 10 plates exhibit better potential to simulate the fresh-snow-layer BRF. In the case of aged snow, shapes of solid and hollow column, hollow bullet rosette, and aggregate of 5 and 10 plates provide BRF values in conformity with respective measurements. However, they slightly underestimate the BRF at high zenith angles in the backscattering region where aggregates of 8 columns can simulate the aged-snow BRF better.

Figure 4 The comparison of angle dependence of laboratory-measured and simulated snow reflectance: (a) fresh snow sample and (b) aged snow sample. Symbols are measurements and lines are simulations with SCIATRAN, assuming different SPSs (see legend).

The above analysis demonstrates that the selected database of SPS can be used successfully to reproduce measured BRF of both fresh and aged snow samples. Similar results were obtained by Pohl et al. (2020b). In this paper, the top-of-atmosphere BRF at 865 nm derived from POLarization and Directionality of the Earth's Reflectances 3 (POLDER-3) on Polarization & Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL) measurements over a pure snow surface in Greenland (70.5^∘ N, 47.3^∘ W) on 6 July 2008 were compared with the SCIATRAN simulations, using droxtals, solid bullet rosettes and solid columns.

According to the above analysis, we can formulate the general algorithm to retrieve SGS and SPS from satellite observations. Satellite provides the wavelength-dependent TOA reflectance. For a given SGS and SPS pair, the minimization between satellite-observed TOA reflectance and theoretical simulation is performed. The optimal SGS and SPS are obtained when the difference between observations and simulations reaches the predefined criteria. The SSA is then calculated by the retrieved SGS and SPS.

4 XBAER algorithm

The retrieval algorithm consists of three stages. The first stage includes the estimation of SGS using the effective Lambertian surface albedo after atmospheric correction for selected observation geometries and wavelengths. This step is performed based on the path radiance representation (Mei et al., 2017), in which the TOA reflectance can be described by the contribution from the atmosphere and the interaction between atmosphere and surface. The inverse to derive the surface reflectance from the satellite-observed TOA reflectance is called the atmospheric correction. And due to certain assumptions in the path radiance representation, the derived surface reflectance is equivalent to the effective Lambertian surface albedo. The estimation of SGS is obtained solving the following minimization problem with respect to the effective radius, r, of snow crystals:

$$\begin{array}{}\text{(1)}& {∥{\mathit{A}}_{\mathrm{e}}-{\mathit{R}}_{\mathrm{s}}\left(r\right)∥}^{\mathrm{2}}\to min.\end{array}$$

Here, A_e and R_s(r) are two vectors whose components are the effective Lambertian surface albedo and the simulated snow reflectance, respectively. The dimension of these vectors is the number of wavelengths times the number of viewing directions.

The simulation of snow reflectance (components of vectors R_s(r)) was performed using the radiative transfer package SCIATRAN (Rozanov et al., 2014) as described in Sect. 3. The optical properties of nine SPSs, listed in Table 1, were used for radiative transfer calculations.

The minimization problem formulated by Eq. (1) was solved separately for each SPS using Brent's method (Brent, 1973). The solution of the minimization problem for each crystal habit is characterized by the following residual:

$$\begin{array}{}\text{(2)}& {\mathrm{\Delta }}_i={∥{\mathit{A}}_{\mathrm{e}}-{\mathit{R}}_{\mathrm{s}}\left(r_i^{\ast }\right)∥}^{\mathrm{2}},i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{9},\end{array}$$

where $r_i^{\ast }$ is the solution of minimization problem given by Eq. (1) for the ith shape of the ice crystal particle.

The second stage is the selection of such i values (SPS) for which Δ_i is minimal. This completes the retrieval process and enables the optimal SGS and SPS to be obtained.

The third stage is to calculate SSA for the retrieved SGS and SPS. To this end, let us rewrite the SSA introduced above in the following equivalent form:

$$\begin{array}{}\text{(3)}& \text{SSA}=\mathrm{3}/\mathit{\rho }r\cdot (A_{\mathrm{t}}/\mathrm{4}{A}_{\mathrm{p}}),\end{array}$$

where r is the effective radius. According to Cauchy's surface area formula (Cauchy, 1841; Tsukerman and Veomett, 2016), the average area of the projections of a convex body is equal to the surface area of the body, up to a multiplicative constant. In our case, this results in A_t=4A_p and SSA for convex particles such as droxtals, solid columns and plates is equal to $\mathrm{3}/\mathit{\rho }r$. In the case of non-convex particles, the calculation of SSA requires information about total area A_t. Although the database given by Yang et al. (2013) does not contain information about A_t, the total area of non-convex particles can be calculated employing the geometric parameters of ice crystal habits presented in Table 1 of Yang et al. (2013). Here we take a typical SPS, aggregate of 8 columns, as an example, to show the difference between SSA calculated assuming convex and non-convex particles.

According to Masanori Saito (personal communication, 2021), the parameters L and a of the aggregate of 8 columns (see Fig. 3 in Yang et al., 2013, for details) can be obtained by scaling with respect to the maximum dimension, D. To find these values for different maximal dimensions, we calculate at first the volume of the aggregate of 8 columns corresponding to parameters a and L on a relative scale as given in Table 1 of Yang et al. (2013).

$$\begin{array}{}\text{(4)}& V_{\mathrm{r}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{\mathrm{8}}{a}_i^{\mathrm{2}}{L}_i\end{array}$$

Using the database of Yang et al. (2013), one can obtain the maximal dimension, D_r, corresponding to the volume, V_r. Introducing the scaling factor, $C_k=D_k/D_{\mathrm{r}}$, we have the semi-width and length for the aggregate with the maximal dimension D_k:

$$\begin{array}{}\text{(5)}& a_{i,k}=a_i{C}_k,L_{i,k}=L_i{C}_k.\end{array}$$

The total surface of the aggregate on a relative scale is given by

$$\begin{array}{}\text{(6)}& S_{\mathrm{r}}=\mathrm{3}\sum _{i=\mathrm{1}}^{\mathrm{8}}\left(\sqrt{\mathrm{3}}{a}_i^{\mathrm{2}}+\mathrm{2}{a}_i{L}_i\right).\end{array}$$

Accounting for Eq. (5), we have

$$\begin{array}{}\text{(7)}& S=C_k^{\mathrm{2}}{S}_{\mathrm{r}}.\end{array}$$

Having obtained the total area, one can calculate SSA as the total surface area of a material per unit of mass:

$$\begin{array}{}\text{(8)}& \text{SSA}={\displaystyle \frac{S}{\mathit{\rho }V}}={\displaystyle \frac{S_{\mathrm{r}}}{\mathit{\rho }{C}_k{V}_{\mathrm{r}}}}.\end{array}$$

Comparing SSA of a convex particle equal to $\mathrm{3}/\mathit{\rho }r$ with the result given by Eq. (8), one can easily notice the difference in SSA calculated from different SPSs using the same SGS. The details of such calculations for other non-convex ice crystal habits are given in the Appendix.

The relationship between SSA and SGS for different SPSs is presented in Fig. 5. According to Fig. 5, an almost inverse linear relationship between SSA and SGS can be found. The lines, representing droxtal, plate and column, overlap, indicating the same SSA for convex particles. For other SPSs with the same SGS, SSA is larger compared to convex faceted particles. SSA is restricted to the range of 0–100 m² kg⁻¹ in this investigation (Picard et al., 2009). For example, for SGS=100 µm, the SSA is 32.7 m² kg⁻¹ for convex faceted particles, whereas SSAs for aggregate of 8 columns, hollow bullet rosette, hollow column, aggregate of 5 plates, aggregate of 10 plates and solid bullet rosette are 44.2, 43.4, 37.7, 74.4, 66.8 and 35.6 m² kg⁻¹, respectively. The relative differences range from 9 %–128 %, depending on the SPS. Taking into account the definition of SSA, one can derive the following relationship between SSA convex and non-convex particles: ${\text{SSA}}_{\text{nc}}={\text{SSA}}_{\mathrm{c}}\cdot (A_{\mathrm{t}}/\mathrm{4}{A}_{\mathrm{p}})$, where subscript c and nc denote convex and non-convex particles, respectively. The obtained results reveal that for all non-convex ice crystals under consideration, $A_{\mathrm{t}}/\mathrm{4}{A}_{\mathrm{p}}>\mathrm{1}$ and the ratio $A_{\mathrm{t}}/\mathrm{4}{A}_{\mathrm{p}}$ weakly depends on the SGS.

Figure 5 Relationship between SGS and SSA for different SPSs. For a better illustration, the realistic range of specific surface area is limited to 100 m² kg⁻¹.

5 Impact of model parameters uncertainty

The accuracy of any retrieval algorithm depends not only on measurement errors but also on the uncertainty in parameters which cannot be retrieved. In our case, such parameters are ice crystal roughness, aerosol contamination and cloud contamination. The impacts of these factors on XBAER-derived SGS and SPS have been investigated and will be discussed in this section. The TOA reflectances at selected channels (0.55 and 1.6 µm) and observation directions for SZA, VZA and RAA of 70, 30 and 135^∘ for nadir and 70, 55 and 135^∘ for oblique, respectively, were calculated using the radiative transfer model SCIATRAN. The details of each scenario will be presented in the corresponding sub-section below.

5.1 Impact of snow particle shape

Since the first stage of the XBAER algorithm is to estimate the SGS assuming a given SPS, it is reasonable to investigate the impact of SPS on the retrieval of SGS. The TOA reflectances of a snow layer at 0.55 and 1.6 µm with the above-given observation geometries were calculated using the following settings for snow layer and atmospheric parameters:

• Snow layer. This consists of ice crystals with SPS set to be severely roughened aggregates of 8 columns and maximal dimensions [100, 300, 500, 700, 1000, 2000, 3000, 5000] µm, which corresponds to SGS [15, 45.1, 75.2, 105.3, 150.4, 300.8, 451.3, 752.1] µm.

• Atmosphere. This is excluded.

The simulated snow reflectances were used as components of vector A_e in Eq. (1). Nine SPSs from the database presented in Yang et al. (2013) are used sequentially in the retrieval process. The atmospheric correction is not performed because the atmosphere is excluded in the forward simulations. This enables avoiding additional errors caused by the atmospheric correction and estimates the pure effect of SPS on the retrieval results. Figure 6 shows the impact of the SPS on SGS retrieval. Different colors and line styles indicate different ice crystals used in the retrieval process. The solid black line represents the retrieved SGS assuming the SPS in the retrieval process is the same as in forward simulations. This line agrees well with the 1:1 line, indicating that the retrieval algorithm has been implemented technically correctly. According to Fig. 6, one can see both underestimation and overestimation of SGS depending on the SPS used in retrieval. However, in most cases, an incorrect SPS leads to an underestimation of SGS. In particular, the maximal effect can be seen when ice crystals of plate shape, rather than of the correct aggregate of 8 columns, is used (solid yellow line). This result can be easily explained coming back to the right panel of Fig. 2. Indeed, one can see that the same reflectance of the snow layer can be obtained using the plate shape, instead an aggregate-of-8-columns shape, with a significantly smaller SGS. These results reveal that the SPS is an important parameter affecting the accuracy of retrieved SGS.

Figure 6 Impact of SPS on the retrieval of SGS.

5.2 Impact of SGS and SPS on SSA

Since the SSA is obtained from the retrieved SGS and SPS, an understanding of how the error in SGS and/or SPS propagates to the SSA will provide helpful information to understand the retrieved SSA. Figure 7 shows the impact of SGS (left) and SPS (right) on XBAER-retrieved SSA. The relative error in SGS, ${\mathit{\epsilon }}_{\mathrm{r}}=(r-r^{\prime })/r$, is propagated to the relative error in SSA as ${\mathit{\epsilon }}_{\text{SSA}}=\mathrm{1}-\mathrm{1}/(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{r}})$, and it is independent of the reference SSA. The left panel of Fig. 7 depicts ε_SSA corresponding to ±0.16 of ε_r. One can see that this results in 19 % and −13.8 % of SSA relative errors, which are presented as the upper and lower error boundaries in the left panel of Fig. 7. The systematical error of ±16 % for SGS was obtained as the maximal relative difference between XBAER-retrieved SGS and both in situ and aircraft-measured SGS (as presented in the companion paper – Mei et al., 2021). This represents the worst case of SGS error propagation into SSA.

Figure 7 Impact of SGS and SPS on the retrieval of SSA. (a) SGS errors – the black line with dots indicates the 0 difference for accurate SGS for aggregate of 8 columns, and the grey area indicates the relative error in SSA introduced by a 16 % error in SGS; (b) SPS selection – different color and line styles indicate different SPSs used in the calculation of SSA while the true SPS is set to be “plate” or other convex particles.

The impact of SPS on SSA is demonstrated in the right panel of Fig. 7. As a reference shape, we have selected in this case the plate, which provides the same SSA as other convex particles. One can see that the SSA of non-convex particles overestimates the SSA of convex particles, which is in line with the results presented in Sect. 4. For instance, for the same SGS, the SSA for aggregate of 8 columns (non-convex particle) is about 3 times larger than that for doxtal (convex particle). Since the assumption of the sphere (convex particle) is used to measure SSA in-field measures (Gallet et al., 2009; Nick Rutter, personal communication. 2021), such as observations from SnowEx, the retrieval results of SSA from XBAER will be systematically larger than field measurements in the case of non-convex particles even if the retrieved and measured SGSs are similar. However, a detailed discussion with respect to uncertainty in the campaign-based measurement is out of the scope of this paper.

5.3 Impact of ice crystal surface roughness

Although the surface roughness of ice crystals is not too severe for snow compared to ice cloud due to basic thermodynamics (Colbeck, 1980, 1983), the ice crystal surface roughness (ICSR), indicating ice crystal surface texture, may still be important for the retrieval of snow properties from optical sensors such as SLSTR. The ICSR has been used as a new variable in model simulation (Järvinen et al., 2018). Retrieval algorithms of ice cloud parameters are frequently based on the assumption that the ice crystal surface is smooth (Kokhanovsky et al., 2019). This assumption can however introduce large uncertainty in the ice cloud retrieval parameters and, as a consequence, lead to misunderstanding the impacts of ice clouds on global climate change (Järvinen et al., 2018). However, this issue has not yet been discussed for snow. In general, ice crystal surfaces are rougher in clouds than in snow layers due to metamorphism processes (Colbeck, 1980, 1983; Ulanowski et al., 2014). The investigation of the impact of ICSR on retrieval of snow properties provides valuable information to understand the XBAER algorithm. The ICSR according to Yang et al. (2013) is defined similarly to the definition suggested by Cox and Munk (1954) for the roughness of the sea surface. A parameter σ describes the degree of ICSR. The σ values 0, 0.03 and 0.5 are for three surface roughness conditions: smooth, moderate roughness and severe roughness. And only the above three values are available in the Yang et al. (2013) database. The snow layer reflectances were used as components of the vector A_e in Eq. (1) in the same way as in Sect. 5.1.

Figure 8 Impact of ice crystal surface roughness (ICSR) on the retrieval of SGS (a), SPS (b) and SSA (c). Different colors indicate different ICSRs used in the retrieval.

Figure 8 shows the impact of ICSR on the retrieved SGS, SPS and SSA. The impact of ICSR on SGS and SSA are relatively small for SGSs smaller than ∼300 µm. Ignoring the impact of roughness leads, in general, to a slight overestimation of SGS and an underestimation of SSA. The absolute errors in SGS and SSA introduced by ICSR range from 0.3 %–3 %, depending on SGS. Due to the inverse almost linear relationship between SSA and SGS, as presented in Fig. 5, for the same SPS, an overestimation of SGS leads to an underestimation of SSA. The slight overestimation can be found if less ICSR is taken into account in retrieval because the snow reflectance with the same SGS and SPS for ICSR=0.5 is larger than for ICSR=0.03 due to the lower asymmetry factor of ice crystals with more roughened surface roughness; thus the same surface reflectance observed by satellite requires a larger SGS for the case with ICSR=0.03 used in retrieval in contrast to ICSR=0.5 used in the forward simulation. However, as can be seen from the right panel of Fig. 8, the XBAER algorithm still retrieves the correct SPS ignoring the impact of roughness.

5.4 Impact of aerosol contamination

The impact of aerosols on the retrieval of snow properties using passive remote sensing can be important because there is limited aerosol information over the cryosphere (Mei et al., 2013a, b, 2020a; Tomis et al., 2015) to perform an accurate atmospheric correction. The use of MERRA-simulated AOT, although of good data quality, will still introduce potential aerosol contamination into the XBAER-derived snow properties. The impact of aerosols on snow property retrieval is much smaller over Arctic regions compared to middle–low latitudes (e.g., Canadian Arctic, Tibetan Plateau) due to large absolute uncertainty in the MERRA-simulated aerosols over middle–low latitudes in wintertime. A detailed comparison of how possible aerosol contamination may affect the retrieved snow properties will be included in the companion paper (Mei et al., 2021). In the companion paper, the comparison between satellite-derived and campaign-measured snow properties all over the world will be included. In order to have a better understanding of aerosol contamination on snow property retrieval, the TOA reflectances were calculated at 0.55 and 1.6 µm with the above-given observation geometries using the following settings:

• Snow layer. This is the same as in Sect. 5.1.

• Atmosphere. Aerosol type is set to be weakly absorbing (Mei et al., 2020b) with AOTs of [0.05, 0.08, 0.11]. Other atmospheric parameters are set according to the Bremen 2D chemical transport model (B2D CTM) for April at 75^∘ N (Sinnhuber et al., 2009). It is worth noticing these three AOT values represent background, average and pollution conditions in the Arctic as suggested by Mei et al. (2020a, b).

Figure 9 Impact of aerosol contamination on the retrieval of SGS (a) SPS (b) and SSA (c). Different colors indicate different AOT used in forward simulations. No atmospheric correction is performed in the retrieval, black dash line is the 1:1 line.

Figure 9 shows the impact of aerosol contamination on the SGS (a), SPS (b) and SSA (c) retrieval. These results are obtained by introducing a 50 % error into AOT at the step of atmospheric correction and can be considered the worst case for impact of aerosol contamination on retrieved SGS, SPS and SSA. The surface reflectances estimated after employing the atmospheric correction were used as components of the vector A_e in Eq. (1). One can see that aerosol introduces systematic underestimation of retrieved SGS for the given scenarios and the magnitude of underestimation increases with the increase in AOT. For a typical background Arctic aerosol condition, with AOT=0.05, aerosol contamination introduces errors in SGS of less than 3 % for SGS≤150 µm and less than 7 % for $\mathrm{150}\le \text{SGS}<\mathrm{300}$ µm. The maximal errors introduced by the aerosol contamination increase to 30 % and 37 % in the case of average and pollution conditions for AOT=0.08 and AOT=0.11, respectively. Please note that the AOT values in the Arctic can be even smaller than 0.05, for instance, AOT over Greenland. Thus, the analysis with respect to aerosol contamination is the worst case for a typical Arctic condition.

Figure 10 Schematic representation of two stages of the retrieval process. (a) Determination of effective radius for each ice crystal form. (b) Selection of optimal SGS–SPS pair.

For the case of AOT=0.05, SPSs have been correctly retrieved for all SGS values, indicating that under a typical Arctic clean condition, the impact of aerosols is not so large as to disturb SPS retrieval. In order to demonstrate the two-stage retrieval process and illustrate the impact of aerosols, let us focus on Fig. 10. To facilitate the presentation, we consider the measurement of reflectance at 1.6 µm for a single observation direction (30^∘) and at 0.55 µm for the difference in reflectance at two observation angles (30 and 55^∘). This enables the avoiding of the minimization process given by Eq. (1) and represents the retrieval process in a simple graphic form. The left panel of Fig. 10 depicts the determination of an effective radius for each ice crystal form, assuming the correct shape is aggregate of 8 columns with an effective radius of 105.4 µm. Solid and dotted lines are the surface reflectance of the snow layer consisting of ice crystals with different forms, and the dashed line is the measured reflectance after the atmospheric correction. The obtained SGSs are in the range 40–120 µm, depending on the selected SPS, and presented in Fig. 10 by solid and dotted vertical lines. In the case of correct SPS selection (aggregate of 8 columns) the retrieved SGS is ∼110 µm. The right panel of Fig. 10 shows the second stage of the retrieval process, namely, the selection of an SPS for which the difference between the measured (dashed line) and simulated (solid black line) value is minimal. In the case under consideration the correct shape is selected with an effective radius of ∼110 µm.

For larger-AOT conditions, an inaccurate selection of SPS occurs for all SGS cases, indicating the remaining aerosol information is strong enough to decouple the aerosol contribution from the snow surface characteristic. Thus, a quality flag of SPS, associated with AOT, should be introduced into the retrieval of real satellite data. It is interesting to see that “solid bullet rosette” is the preferable SPS for very strong aerosol contamination cases. This is due to the similar scattering properties (shape) of ice crystals and weakly absorbing aerosols, defined in forward simulation. The impact of aerosol contamination, for typical Arctic conditions, introduces a less-than-5 % error into SSA. However, for large aerosol contamination, the around 30 % underestimation in SGS linearly introduced about 25 % overestimation in SSA, which agrees with the analysis as presented in Fig. 7.

6 Impact of cloud contamination

Any cloud-screening method, especially over the cryosphere, may introduce cloud contamination for the retrieval of atmospheric and surface properties (Chen et al., 2014; Mei et al., 2017; Jafariserajehlou et al., 2019). Understanding of the cloud contamination will provide valuable information to interpret the retrieval results using the SLSTR instrument. To investigate the impact of cloud contamination, the following settings were used to perform the simulations of TOA reflectance:

• Snow layer. This is the same as in Sect. 5.1.

• Atmosphere. An aerosol-free atmosphere is used with other parameters as in Sect. 5.4. Additionally, vertically homogeneous ice cloud consisting of aggregates of 8 columns with an effective radius of 45 µm and optical thickness [0.1, 0.5, 1.0, 5] is set to be at the position of [5, 6] km.

Figure 11 Impact of cloud contamination on the retrieval of SGS (a), SPS (b) and SSA (c). Different colors indicate different COTs in forward simulations; dashed black line is the 1:1 line.

Figure 11 shows the impact of cloud contamination on XBAER-retrieved SGS (a), SPS (b) and SSA (c). The size of ice crystals in ice clouds is typically smaller than snow grain size (Kikuchi et al., 2013). Our statistical analysis of the ice crystal effective radius over Greenland shows an average value in the range of 30–50 µm, which is consistent with previous publications (King et al., 2013; Platnick et al., 2017). According to Fig. 11, an overestimation of SGS can be found for SGSs of less than 45 µm (cloud effective radius) and an underestimation of SGS can be found for SGSs larger than 45 µm. The magnitude of overestimation/underestimation increases with the increase in cloud optical thickness (COT). XBAER-derived SGS becomes saturated for COT larger than 0.5. Due to a limited photon penetration depth for optically thicker clouds (e.g., COT=5), the XBAER algorithm retrieves the effective radius of ice crystals in the cloud. This demonstrates that theoretically, the XBAER algorithm can retrieve an ice cloud effective radius without pre-processing of cloud screening. And this can be further used as post-processing to avoid cloud contamination.

The impact of the cloud on the retrieval of SPS is similar to the impact of aerosols considered above. In short, the cloud plays a larger role for larger SPSs (darker TOA) and this impact increases with the increase in COT. However, cloud with large COT can be much more easily detected and excluded by the cloud-screening algorithm (e.g., for the cases with COT>0.5). SPSs are correctly picked up due to the same SPS used for both the snow layer and the cloud layer. Similarly to the impact of aerosol, the underestimation of SGS introduced by the cloud leads to an overestimation of SSA (Fig. 11c). The increase in COT results in saturation of the ice cloud SSA, with a value of 100 m² kg⁻¹ in the case of aggregates of 8 columns.

7 Impact of other factors occurring in reality

The above theoretical investigations include all possible important factors affecting the accuracy of the XBAER algorithm. However, when applying the XBAER algorithm to the SLSTR instrument for real scenarios, two additional factors need to be considered as well: one is the impact of the instrument spectral response function (SRF), and the other one is the representativeness of the snow scenario for reality.

7.1 Impact of instrument spectral response function

• Snow layer. This is the same as in Sect. 5.1.

• Atmosphere. An aerosol-free atmosphere is used with other parameters as in Sect. 5.4.

The forward simulations are performed with and without the impact of the spectral response function (SRF). The SRFs for SLSTR at 0.55 and 1.6 µm are shown in Fig. 12. The retrieval is then performed ignoring the SRF. Figure 13 shows the impact of SRF on the retrieval of SGS, SPS, and SSA. For forward simulations without taking the SRF into account (labeled as No in Fig. 13), SGS, SPS and SSA are well received as expected. And this agrees with Fig. 6. However, ignoring the impact of the SRF introduces about 7 % uncertainties into the simulated surface reflectance, and this causes about a 5 %–7 % error in both SGS (overestimation) and SSA (underestimation). Taking the SRF into account leads to a smaller surface reflectance at 1.6 µm due to potential gas absorption at this wavelength and thus introduces an overestimation for SGS. However, due to a significantly smaller impact at 0.55 µm, the SRF does not play a significant role in the retrieval of SPS.

Figure 12 Spectral response function of (a)  0.55 and (b) 1.6 µm of the SLSTR instrument.

Figure 13 Impact of SRF on the retrieval of SGS (a), SPS (b) and SSA (c). Different colors indicate retrieval results without (No) and with (Yes) the SRF in forward simulations; dashed black line is the 1:1 line.

7.2 Impact of snow inhomogeneities

In this section, a realistic model of the snow layer is represented by a vertically inhomogeneous, polydisperse ice crystal habit mixture. Following Saito et al. (2019), the gamma distribution with respect to the maximal dimension will be used to describe polydisperse properties:

$$\begin{array}{}\text{(9)}& n\left(D\right)=NG\left(D\right).\end{array}$$

Here, N is the number of ice particles per unit volume and G(D) is the gamma distribution function; i.e.,

$$\begin{array}{}\text{(10)}& G\left(D\right)=C(D/v{)}^{k-\mathrm{1}}{e}^{-D/v},\end{array}$$

where k and v are the shape and scale parameters, the normalization factor C is defined as

$$\begin{array}{}\text{(11)}& C={\left[\underset{D_{min}}{\overset{D_{max}}{\int }}(D/v{)}^{k-\mathrm{1}}{e}^{-D/v}\mathrm{d}D\right]}^{-\mathrm{1}},\end{array}$$

and D_min and D_max describe the minimal and maximum particle sizes in the distribution.

In order to introduce the vertical inhomogeneity, we use the measurement of snow density and equivalent optical diameter vertical profiles conducted during the SnowEx17 campaign. Accounting for the fact that the equivalent optical diameter cannot be directly used to define parameters of the gamma distribution, we use the vertical profile as a shape of the mode (most frequent value in a dataset); i.e.,

$$\begin{array}{}\text{(12)}& D_{\mathrm{0}}\left(z\right)={\displaystyle \frac{D_{\mathrm{e}}\left(z\right)}{D_{\mathrm{e}}\left(z_{\text{top}}\right)}}{D}_{\mathrm{0}}\left(z_{\text{top}}\right),\end{array}$$

where D_e(z) is the measured vertical profile of equivalent optical diameter and D₀(z) is the vertical profile of the mode. The mode near the top of the snow layer, D₀(z_top), is assumed to be equal to 400 µm according to the measurement data reproduced by Saito et al. (2019) in their Fig. A1.

Taking into account the analytical expression of the mode via shape and scale parameters,

$$\begin{array}{}\text{(13)}& D_{\mathrm{0}}=(k-\mathrm{1})v,\end{array}$$

and the relationship between shape and scale parameters derived by Saito et al. (2019),

$$\begin{array}{}\text{(14)}& k=\mathrm{11.38}{v}^{-\mathrm{0.167}}-\mathrm{2},\end{array}$$

we can estimate parameters k and v of the gamma distribution corresponding to D₀(z) given by Eq. (12).

The snow grain habit mixture (SGHM) model is used according to Saito et al. (2019). In particular, the particle habits include droxtal, solid hexagonal column and solid bullet rosette. The habit fraction, f_h(D), as a function of the maximal dimension of the SGHM model is presented in Fig. 15b. The habit fraction is defined so that, for each D,

$$\begin{array}{}\text{(15)}& \sum _{h=\mathrm{1}}^{\mathrm{3}}{f}_{\mathrm{h}}\left(D\right)=\mathrm{1}.\end{array}$$

The selected SGHM model enables us to derive the total volume of ice per unit volume of air as

$$\begin{array}{}\text{(16)}& V_{\mathrm{t}}=N\sum _{h=\mathrm{1}}^{\mathrm{3}}\left[\underset{D_{min}}{\overset{D_{max}}{\int }}{V}_h\left(D\right)f_{\mathrm{h}}\left(D\right)G\left(D\right)\mathrm{d}D\right]\end{array}$$

and ice water content (IWC)

$$\begin{array}{}\text{(17)}& \text{IWC}=V_{\mathrm{t}}{\mathit{\rho }}_{\text{ice}},\end{array}$$

where V_h(D) is the volume of each habit as given in the database of Yang et al. (2013) and ρ_ice is the density of ice.

Figure 14 Snow properties used for simulations to investigate the impacts of snow layer model on the XBAER retrieval (a) snow grain size profile and (b) snow density observed during the SnowEx17 campaign.

Taking into account the vertical profile if IWC is measured (see right panel of Fig. 14), we can obtain the vertical profile of the particle number density. Using Eqs. (16) and (17), we have

$$\begin{array}{}\text{(18)}& N\left(z\right)={\displaystyle \frac{\text{IWC}\left(z\right)}{{\mathit{\rho }}_{\text{ice}}\sum _{h=\mathrm{1}}^{\mathrm{3}}\left[\underset{D_{min}}{\overset{D_{max}}{\int }}{V}_h\left(D\right)f_{\mathrm{h}}\left(D\right)G(D,z)\mathrm{d}D\right]}}.\end{array}$$

Summing up, we define the microphysical properties of the snow layer using the following model of particle size distribution:

$$\begin{array}{}\text{(19)}& n(D,z)=N\left(z\right)C{\left[{\displaystyle \frac{(\stackrel{\mathrm{‾}}{k}-\mathrm{1})D}{D_{\mathrm{0}}\left(z\right)}}\right]}^{\stackrel{\mathrm{‾}}{k}-\mathrm{1}}\exp \left[-{\displaystyle \frac{(\stackrel{\mathrm{‾}}{k}-\mathrm{1})D}{D_{\mathrm{0}}\left(z\right)}}\right],\end{array}$$

where D₀(z) and N(z) are given by Eqs. (12) and (18), respectively, and the shape parameter, k, is assumed to be altitude independent and set to 2.3.

Figure 15 Snow properties used for simulations to investigate the impacts of the habit mixture model on XBAER retrieval: (a) particle size distribution of snow grain size in snow layer and (b) habit fraction suggested by Saito et al. (2019).

The bulk single-scattering properties of the snow layer such as the extinction coefficient, scattering coefficient and scattering function are defined in the same way as proposed by Baum et al. (2011). For instance, the bulk extinction coefficient is calculated as

$$\begin{array}{}\text{(20)}& {\mathit{\beta }}_{\text{ext}}\left(z\right)=\underset{D_{min}}{\overset{D_{max}}{\int }}\left[\sum _{h=\mathrm{1}}^{\mathrm{3}}{\mathit{\sigma }}_{\text{ext},h}\left(D\right)f_{\mathrm{h}}\left(D\right)n(D,z)\mathrm{d}D\right],\end{array}$$

where σ_ext,h(D) is the extinction cross section as given for each habit in the database of Yang et al. (2013).

The following settings are used to simulate the reference snow reflectance at wavelengths 0.55 and 1.6 µm:

• Snow layer. This is a vertically inhomogeneous, polydisperse habit mixture and model as described above.

• Atmosphere. This is excluded.

Using the simulated reflectances in the XBAER algorithm, we have retrieved SPS as droxtal with the maximal dimension equal to 740 µm. Taking into account that the model of particle size distribution (PSD) near the top of the snow layer is 400 µm and the mean value calculated as $k{D}_{\mathrm{0}}/(k-\mathrm{1})$ is equal to 708 µm, one can see that the retrieved maximal dimension is an estimation of the mean value of PSD near the top of the snow layer.

Since there is no single reference SGS value when a PSD is used, it is important to check the representativeness of XBAER-derived SGS. Accounting for the fact that the mode and mean values for a given PSD are two typical “effective” ways to describe a polydisperse medium, we compared reflectances of the snow layer calculated assuming PSD in the form of a gamma distribution and assuming a monodisperse medium with SGS equal to the mode or to the mean of the selected PSD. In order to simplify analysis, we consider vertically homogeneous snow layer consisting of only single-particle habit. The calculations of reflectance were performed for severely roughened aggregate-of-8-column and droxtal particles setting the shape parameter, k, equal to 2.3 and the model equal to [100, 300, 500, 700, 1000, 2000, 3000, 5000] µm.

Figure 16 The comparison between simulated snow reflectance using a mono- and polydisperse snow model consisting of aggregate of columns (a) and droxtal (b). In the case of the monodisperse model, the SGS is assumed to be the mean and mode value of PSD at the top of the snow layer (see Fig. 14a). The reference value is shown on the x axis.

Figure 16 shows the comparison of snow reflectance calculated assuming a monodisperse and polydisperse snow model. In the case of the monodisperse model, SGS is assumed to be equal to the mean or to the mode value of PSD. We can see that the surface reflectance calculated using the mean value of PSD agrees better with reference values than reflectance calculated using the mode value. In particular, the root-mean-square error (RMSE) values are more than 2 times smaller. One can also see from Fig. 16 that the difference between monodisperse reflectances calculated using mean or mode PSD values decreases with an increase in the PSD mode. This can be explained due to the fact that the increase in the PSD mode leads to the increase in absorption and decrease in reflectance sensitivity with respect to the variation in SGS.

8 Conclusions

SGS, SPS and SSA are three important parameters to describe snow properties. They play important roles in the changes in snow albedo/reflectance and impact the atmospheric and energy exchange processes. A better knowledge of SGS, SPS and SSA can provide more accurate information to describe the impact of snow on Arctic amplification processes. The information about SGS, SPS and SSA may also be used explore new applications to understand atmospheric conditions (e.g., aerosol loading). Although some previous attempts (e.g., Lyapustin et al., 2009) show the capabilities of using passive remote sensing to derive SGS over a large scale, no publication has been found in which SGS, SPS and SSA are derived simultaneously. To our best knowledge, this is the first paper attempting to retrieve these parameters using satellite observations.

The new algorithm is designed within the framework of the XBAER algorithm. The XBAER algorithm has been applied to derive SGS, SPS and SSA using the newly launched SLSTR instrument on board the Sentinel-3 satellite. This is the first part of two companion papers to describe the algorithm and to present the sensitivity studies.

The SGS, SPS and SSA retrieval algorithm is based on the recent publication by Yang et al. (2013), in which a database of optical properties for nine typical SPSs (aggregate of 8 columns, droxtal, hollow bullet rosette, hollow column, plate, aggregate of 5 plates, aggregate of 10 plates, solid bullet rosette, column) are provided. Previous publications show that this database can be used to retrieve ice crystal properties in both ice cloud and snow (e.g., Järvinen et al., 2018; Saito et al., 2019). The algorithm is an LUT-based approach, in which the minimization is achieved by the comparison between atmospherically corrected TOA reflectance at 0.55 and 1.6 µm observed by SLSTR and a pre-calculated LUT of surface reflectances under different geometries and snow properties. The retrieval is relatively time-consuming because the minimization has to be performed for each SPS and the optimal SGS and SPS are selected after nine minimizations are performed. The SSA is then obtained using the retrieved SGS and SPS based on another pre-calculated LUT.

The sensitivity studies with respect to the impacts of SPS, ICSR, and aerosol and cloud contamination on XBAER-derived SGS and SPS provide a comprehensive understanding of the retrieval accuracy of the new algorithm. The main findings of the theoretical considerations are (1) XBAER-derived SGS is more likely to represent the average SGS near the top of snow layer when a PSD is known; (2) SPS plays an important role in the retrieval accuracy of SGS – the retrieved SGS can differ several times by usage different SPSs in the retrieval process; (3) impact of ICSR on the retrieval accuracy of SGS can be neglected and ignoring ICSR completely may introduce a maximal 3 % error into the retrieval accuracy of SGS, especially for large ice crystals; (4) the assumption of a convex particle shape (e.g., sphere) for a non-convex ice crystal leads to the underestimation of the retrieved SSA; (5) the impact of aerosols and cloud increases with the increase in both aerosol/cloud loading and SGS; and (6) the impact of instrument SRF may introduce some positive bias for SGS and negative bias for SSA; however, it plays no role for the determination of SPS.

Even though all major possible factors affecting the retrieval accuracy of the XBAER algorithm are investigated in this paper, in reality, the final retrieval accuracy can only be evaluated by performing a thorough comparison with independent measurement results because uncertainties caused by each individual factor can compensate for each other in the real satellite retrieval. All details of such validation can be found in the companion paper of Mei et al. (2021).

Appendix A

According to the definition of specific surface area,

$$\begin{array}{}\text{(A1)}& \text{SSA}={\displaystyle \frac{A}{\mathit{\rho }V}},\end{array}$$

one needs to calculate the total area A of ice crystals. In the following sections, we consider in details the basic equations to calculate total area and SSA of different SPSs given in the database of Yang et al. (2013) and used above within the retrieval algorithm.

A1 Droxtal, solid column, plate

In the case of convex faceted particles such as droxtal, solid column and plate, the calculation of total area is straightforward and based on Cauchy's surface area formula:

$$\begin{array}{}\text{(A2)}& A=\mathrm{4}{A}_{\mathrm{p}}.\end{array}$$

Taking into account that for selected SPSs, one can find corresponding V and A_p values in the database given by Yang et al. (2013), we have the following result for the SSA of such particles:

$$\begin{array}{}\text{(A3)}& \text{SSA}={\displaystyle \frac{\mathrm{4}{A}_{\mathrm{p}}}{\mathit{\rho }V}}.\end{array}$$

A2 Hollow column

In this case a solid column includes two equal cavities in the form of a hexagonal pyramid and cannot be considered a convex particle. The aspect ratio of the hollow column with the height, d, of a hexagonal pyramid is given according to Yang et al. (2013) as

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{2}a}{L}}=\left\{\begin{array}{ll}\mathrm{0.7},& L<\mathrm{100}\mathrm{\mu }\mathrm{m}\\ \frac{\mathrm{6.96}}{\sqrt{L}},& L\ge \mathrm{100}\mathrm{\mu }\mathrm{m}\end{array}\right.,d=\mathrm{0.25}L.\end{array}$$

The volume of such a hollow column is given by

$$\begin{array}{}\text{(A5)}& V=V_{\mathrm{c}}-\mathrm{2}{V}_{\mathrm{p}},\end{array}$$

where the volume of the solid column, V_c, and a hexagonal pyramid, V_p, are

$$\begin{array}{}\text{(A6)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{c}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}L,\text{(A7)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{p}}={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}d.\end{array}$$

Thus, the volume, V, is

$$\begin{array}{}\text{(A8)}& V={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}(\mathrm{3}L-\mathrm{2}d).\end{array}$$

Employing the relationship between d and L given by Eq. (A4) and excluding a, we have

$$\begin{array}{}\text{(A9)}& V={\displaystyle \frac{\mathrm{2.5}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}L=\left\{\begin{array}{ll}{m}_{\mathrm{0}}{m}_{\mathrm{1}}^{\mathrm{2}}{L}^{\mathrm{3}},& L<\mathrm{100}\mathrm{\mu }\mathrm{m}\\ m_{\mathrm{0}}{m}_{\mathrm{2}}^{\mathrm{2}}{L}^{\mathrm{2}},& L\ge \mathrm{100}\mathrm{\mu }\mathrm{m}\end{array}\right.,\end{array}$$

where m₀ is $\mathrm{2.5}\times \sqrt{\mathrm{3}}/\mathrm{2}$, m₁ is $\frac{\mathrm{0.7}}{\mathrm{2}}$ and m₂ is $\frac{\mathrm{6.96}}{\mathrm{2}}$. For a selected volume, V, the length, L, is calculated as follows:

$$\begin{array}{}\text{(A10)}& L=\left\{\begin{array}{ll}{\left[V/m_{\mathrm{0}}/m_{\mathrm{1}}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}},& V<V_{\mathrm{100}}\\ {\left[V/m_{\mathrm{0}}/m_{\mathrm{2}}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},& V\ge V_{\mathrm{100}}\end{array}\right.,\end{array}$$

where V₁₀₀ is $m_{\mathrm{0}}{m}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{100}}^{\mathrm{2}}.$

Let us now calculate the area of each triangle side of the pyramid:

$$\begin{array}{}\text{(A11)}& S_{\mathrm{t}}={\displaystyle \frac{a}{\mathrm{2}}}\sqrt{d^{\mathrm{2}}+{\displaystyle \frac{\mathrm{3}{a}^{\mathrm{2}}}{\mathrm{4}}}}.\end{array}$$

The area of the lateral surface of two pyramids is

$$\begin{array}{}\text{(A12)}& S_{\mathrm{p}}=\mathrm{3}a\sqrt{\mathrm{4}{d}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}.\end{array}$$

And the total surface area of the hollow column is given by

$$\begin{array}{}\text{(A13)}& S=\mathrm{6}aL+\mathrm{3}a\sqrt{\mathrm{4}{d}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}},\end{array}$$

where a and d should be expressed via L according to Eq. (A4).

Having obtained the total area, one can calculate the specific surface area:

$$\begin{array}{}\text{(A14)}& \text{SSA}={\displaystyle \frac{S}{\mathit{\rho }V}}.\end{array}$$

A3 Hollow bullet rosettes

In this case a solid column includes a cavity in the form of a hexagonal pyramid with height H and a hexagonal pyramid with height t on the opposite site of the column. The aspect ratio and parameters H and t are given according to Yang et al. (2013) as

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{2}a}{L}}=& \mathrm{2.3104}{L}^{-\mathrm{0.37}},t={\displaystyle \frac{\sqrt{\mathrm{3}}a}{\mathrm{2}\tan \left({\mathrm{28}}^{{}^{\circ }}\right)}},\\ & H=\mathrm{0.5}(t+L).\end{array}\end{array}$$

The volume of a hollow bullet rosette is given by

$$\begin{array}{}\text{(A16)}& V_{\mathrm{1}}=V_{\mathrm{c}}-V_{-}+V_{+}.\end{array}$$

Using Eqs. (A6) and (A7), we have

$$\begin{array}{ll}\text{(A17)}& {\displaystyle }{V}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}L-{\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}H+{\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}t{\displaystyle }& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}(\mathrm{3}L-H+t).\end{array}$$

Substituting H as given by Eq. (A15), we obtain

$$\begin{array}{}\text{(A18)}& V_{\mathrm{1}}={\displaystyle \frac{\sqrt{\mathrm{3}}{a}^{\mathrm{2}}}{\mathrm{4}}}(\mathrm{5}L+t).\end{array}$$

Using Eq. (A15), we express parameters a and t of hollow bullet rosettes via L:

$$\begin{array}{}\text{(A19)}& {\displaystyle }& {\displaystyle }a=m_{\mathrm{a}}{L}^{p_{\mathrm{a}}},\text{(A20)}& {\displaystyle }& {\displaystyle }t=m_{\mathrm{t}}{m}_{\mathrm{a}}{L}^{p_{\mathrm{a}}},\end{array}$$

where coefficients m_a, m_t and p_a are

$$\begin{array}{}\text{(A21)}& m_{\mathrm{a}}={\displaystyle \frac{\mathrm{2.3104}}{\mathrm{2}}},m_{\mathrm{t}}={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}\tan \left({\mathrm{28}}^{{}^{\circ }}\right)}},p_{\mathrm{a}}=\mathrm{1}-\mathrm{0.37}.\end{array}$$

The expression Eq. (A18) can be rewritten as

$$\begin{array}{}\text{(A22)}& V_{\mathrm{1}}={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{4}}}{m}_{\mathrm{a}}^{\mathrm{2}}{L}^{\mathrm{2}{p}_{\mathrm{a}}+\mathrm{1}}(\mathrm{5}+m_{\mathrm{t}}{m}_{\mathrm{a}}{L}^{-\mathrm{0.37}}).\end{array}$$

The total area of a hollow bullet rosette is written as

$$\begin{array}{}\text{(A23)}& S_{\mathrm{1}}=\mathrm{6}aL+{\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}\sqrt{\mathrm{4}{H}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}\sqrt{\mathrm{4}{t}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}\end{array}$$

and can be calculated when for a selected maximal dimension D the parameter L is known. For a desired dimension D (volume V) of hollow bullet rosettes, consisting of six equal rosettes (see Table 1), Eq. (A22) was solved with respect to the length, L, using the following iterative approach:

$$\begin{array}{}\text{(A24)}& L_n={\left[{\displaystyle \frac{\mathrm{2}V}{\mathrm{3}\sqrt{\mathrm{3}}{m}_{\mathrm{a}}^{\mathrm{2}}(\mathrm{5}+m_{\mathrm{t}}{m}_{\mathrm{a}}{L}_{n-\mathrm{1}}^{-\mathrm{0.37}})}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}{p}_{\mathrm{a}}+\mathrm{1}}}.\end{array}$$

The iterative process starts with L₀=1 and finishes when $\left|\right(L_n-L_{n-\mathrm{1}})/L_n|\le {\mathrm{10}}^{-\mathrm{4}}$. The total area of hollow bullet rosettes is calculated as

$$\begin{array}{}\text{(A25)}& S_{\mathrm{1}}=\mathrm{6}aL+{\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}\sqrt{\mathrm{4}{H}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}\sqrt{\mathrm{4}{t}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}.\end{array}$$

The SSA is given by

$$\begin{array}{}\text{(A26)}& \text{SSA}={\displaystyle \frac{\mathrm{6}{S}_{\mathrm{1}}}{\mathit{\rho }V}}.\end{array}$$

A4 Solid bullet rosettes

The aspect ratio and parameter t are given according to Yang et al. (2013) as

$$\begin{array}{ll}\text{(A27)}& {\displaystyle }& {\displaystyle \frac{\mathrm{2}a}{L}}=\mathrm{2.3104}{L}^{-\mathrm{0.37}},{\displaystyle }& {\displaystyle }t={\displaystyle \frac{\sqrt{\mathrm{3}}a}{\mathrm{2}\tan \left({\mathrm{28}}^{{}^{\circ }}\right)}}.\end{array}$$

The volume of a single solid bullet rosette is

$$\begin{array}{}\text{(A28)}& V_{\mathrm{1}}=V_{\mathrm{c}}+V_{+}.\end{array}$$

Using Eq. (A6), we have

$$\begin{array}{ll}\text{(A29)}& {\displaystyle }{V}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}L+{\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}t{\displaystyle }& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}^{\mathrm{2}}(\mathrm{3}L+t).\end{array}$$

Using the formula given by Eq. (A27), we express parameters a and t of a solid bullet rosette via L:

$$\begin{array}{}\text{(A30)}& {\displaystyle }& {\displaystyle }a=m_{\mathrm{a}}{L}^{p_{\mathrm{a}}},\text{(A31)}& {\displaystyle }& {\displaystyle }t=m_{\mathrm{t}}{m}_{\mathrm{a}}{L}^{p_{\mathrm{a}}},\end{array}$$

where coefficients, m_a, m_t and p_a are the same as in the case of hollow bullet rosettes given by Eq. (A21). The expression Eq. (A29) can be rewritten as

$$\begin{array}{}\text{(A32)}& V_{\mathrm{1}}={\displaystyle \frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}{m}_{\mathrm{a}}^{\mathrm{2}}{L}^{\mathrm{2}{p}_{\mathrm{a}}+\mathrm{1}}(\mathrm{3}+m_{\mathrm{t}}{m}_{\mathrm{a}}{L}^{-\mathrm{0.37}}).\end{array}$$

For a desired volume V of solid bullet rosettes, consisting of six equal rosettes (see Table 1), this equation was solved with respect to the length, L, of the solid bullet rosette using the following iterative approach:

$$\begin{array}{}\text{(A33)}& L_n={\left[{\displaystyle \frac{V}{\mathrm{3}\sqrt{\mathrm{3}}{m}_{\mathrm{a}}^{\mathrm{2}}(\mathrm{3}+m_{\mathrm{t}}{m}_{\mathrm{a}}{L}_{n-\mathrm{1}}^{-\mathrm{0.37}})}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}{p}_{\mathrm{a}}+\mathrm{1}}}.\end{array}$$

The total area of solid bullet rosettes is calculated as

$$\begin{array}{}\text{(A34)}& S_{\mathrm{1}}=\mathrm{6}aL+{\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}{a}^{\mathrm{2}}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{3}a}{\mathrm{2}}}\sqrt{\mathrm{4}{t}^{\mathrm{2}}+\mathrm{3}{a}^{\mathrm{2}}}.\end{array}$$

The SSA is given by

$$\begin{array}{}\text{(A35)}& \text{SSA}={\displaystyle \frac{\mathrm{6}{S}_{\mathrm{1}}}{\mathit{\rho }V}}.\end{array}$$

A5 Aggregate of 5 and 10 plates

In accordance with the paper of Yang et al. (2013), Table 1 provides the aspect ratios of the ice crystal habits. In the case of an aggregate of columns or plates, the semi-width a and length L of each hexagonal element of the aggregate are on a relative scale. In order to convert these parameters into absolute values, let us consider the following relationship given in Yang et al. (2013) for the aspect ratio of the plate:

$$\begin{array}{}\text{(A36)}& {\displaystyle \frac{\mathrm{2}a}{L}}=\left\{\begin{array}{ll}\mathrm{1},& a\le \mathrm{2}\mathrm{\mu }\mathrm{m}\\ m_{\mathrm{1}}a+m_{\mathrm{0}},& \mathrm{2}<a<\mathrm{5}\mathrm{\mu }\mathrm{m}\\ m{a}^p,& a\ge \mathrm{5}\mathrm{\mu }\mathrm{m}\end{array}\right.,\end{array}$$

where constants are m₁=0.2914, m₀=0.4172, m=0.8038 and p=0.526.

Using this expression and accounting for the fact that relative values for a, given in Table 1, are greater than 5 µm, we can express L_r via a_r as

$$\begin{array}{}\text{(A37)}& L_{\mathrm{r}}={\displaystyle \frac{\mathrm{2}{a}_{\mathrm{r}}}{m{a}_{\mathrm{r}}^p}}={\displaystyle \frac{\mathrm{2}{a}_{\mathrm{r}}^{\mathrm{1}-p}}{m}},\end{array}$$

where subscript r denotes that they are on relative scale. The volume of a hexagonal plate on a relative scale is given by

$$\begin{array}{}\text{(A38)}& v_{\mathrm{r}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}_{\mathrm{r}}^{\mathrm{2}}{L}_{\mathrm{r}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{m}}{a}_{\mathrm{r}}^{\mathrm{3}-p}.\end{array}$$

The volume of aggregates of 5 or 10 plates is given by

$$\begin{array}{}\text{(A39)}& V_{\mathrm{r}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{m}}\sum _{i=\mathrm{1}}^N{a}_{\mathrm{r},i}^{\mathrm{3}-p},\end{array}$$

where N is 5 and N is 10 for 5 and 10 plates, respectively. The absolute value of the volume, V, for a selected maximal dimension of aggregates of 5 or 10 plates can be found in the database presented by Yang et al. (2013). Introducing the scaling factor,

$$\begin{array}{}\text{(A40)}& C={\displaystyle \frac{V_{\mathrm{r}}}{V}},\end{array}$$

we rewrite expression Eq. (A38) as

$$\begin{array}{}\text{(A41)}& V={\displaystyle \frac{V_{\mathrm{r}}}{C}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{mC}}\sum _{i=\mathrm{1}}^N{a}_{\mathrm{r},i}^{\mathrm{3}-p}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{m}}\sum _{i=\mathrm{1}}^N{a}_i^{\mathrm{3}-p},\end{array}$$

where the absolute value of semi-width, a_i, is given by

$$\begin{array}{}\text{(A42)}& a_i={\displaystyle \frac{a_{\mathrm{r},i}}{C^{(\mathrm{3}-p{)}^{-\mathrm{1}}}}}.\end{array}$$

Having obtained the absolute value of a_i for each plate, the length L_i is calculated as

$$\begin{array}{}\text{(A43)}& L_i=\left\{\begin{array}{ll}\mathrm{2}{a}_i,& a\le \mathrm{2}\mathrm{\mu }\mathrm{m}\\ \frac{\mathrm{2}{a}_i}{m_{\mathrm{1}}{a}_i+m_{\mathrm{0}}},& \mathrm{2}<a<\mathrm{5}\mathrm{\mu }\mathrm{m}\\ \frac{\mathrm{2}{a}^{(\mathrm{1}-p)}}{m},& a\ge \mathrm{5}\mathrm{\mu }\mathrm{m}\end{array}\right..\end{array}$$

The total area of a hexagonal plate with semi-width a_i and length L_iis given by

$$\begin{array}{}\text{(A44)}& S_i=\mathrm{2}{\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}}{a}_i^{\mathrm{2}}+\mathrm{6}{a}_i{L}_i.\end{array}$$

The total area is given by

$$\begin{array}{}\text{(A45)}& S=\sum _{i=\mathrm{1}}^N{S}_i.\end{array}$$

Having obtained the total area, one can calculate SSA as the total surface area of a material per unit of mass:

$$\begin{array}{}\text{(A46)}& \text{SSA}={\displaystyle \frac{S}{\mathit{\rho }V}},\end{array}$$

where ρ=917 kg m⁻³ is the density of ice.