Parameterization of single-scattering properties of snow

Introduction Conclusions References Tables Figures


Introduction
Snow grains are non-spherical and often irregular in shape. Still, in many studies, spherical snow grains have been assumed in radiative transfer calculations due to the convenience of using Mie theory. In fact, it has been shown that the spectral albedo of snow can be fitted by radiative transfer calculations under the assumption of spherical snow grains, when the effective snow grain size is considered an adjustable parameter (i.e. determined based on albedo rather than microphysical measurements) (Wiscombe and Warren, 1980;Grenfell et al., 1994;Aoki et al., 2000). Snow albedo parameterizations used in climate models and numerical weather prediction models are often semi-empirical and do not specify the snow grain shape (for some examples, see Wang and Zeng, 2010). However, in most (if not all) physically based albedo parameterizations that explicitly link the albedo to snow grain size, spherical snow grains are assumed (Flanner and Zender, 2005;Gardner and Sharp, 2010;Aoki et al., 2011).
It is, however, well known that the single-scattering properties (SSPs) of non-spherical particles, including the singlescattering albedo ω, the phase function P 11 and the entire phase matrix P, can differ greatly from those of spheres. 1 A consequence of this is that the assumed shape of snow grains has a profound effect on the bidirectional reflectance distribution function (BRDF) of snow (Mischenko et al., 1999;Xie et al., 2006). Furthermore, Aoki et al. (2000) showed that the modelled BRDF of snow agreed better with observa-P. Räisänen et al.: Parameterization of single-scattering properties of snow tions when, instead of the actual phase function for spheres, the Henyey-Greenstein (HG) phase function (Henyey and Greenstein, 1941) was assumed. The HG phase function is very smooth, while that of spheres features icebow and glory peaks not seen for real snow along with very low sideward scattering. Based on a comparison of a few shape models with phase function measurements for laboratory-generated ice crystals (Barkey et al., 2002), Kokhanovsky and Zege (2004) recommended, instead of spheres, the use of Gaussian random spheres (Muinonen et al., 1996;Nousiainen and Muinonen, 1999) or Koch fractals (Macke et al., 1996), which both exhibit a relatively featureless phase function. Since Gaussian random spheres have several free parameters while Koch fractals have none (except for the degree of distortion for randomized Koch fractals), Koch fractals were selected by Kokhanovsky and Zege (2004). Kokhanovsky et al. ( , 2011 further demonstrated that the reflectance patterns computed for Koch fractals agreed reasonably well with actual measurements for snow. Subsequently, they have been used in several studies related to remote sensing of snow grain size and snow albedo (Lyapustin et al., 2009;Negi and Kokhanovsky, 2011;Kokhanovsky et al., 2011).
Other snow grain shape models have also been considered. Tanikawa et al. (2006) suggested the use of non-spherical ice particles with rough surfaces, specifically, cylindrical particles for new snow and prolate ellipsoids for old granular snow. These choices improved the agreement with observed angular reflectance patterns compared to the use of spheres. Jin et al. (2008) compared anisotropic reflectance factors computed using spheres, hexagonal plates, hexagonal columns and aggregates of columns with ground-based measurements in Antarctica, finding the best agreement for the aggregate model and the largest discrepancies for spheres. Furthermore, Zege et al. (2011) tested, in their retrieval algorithm of snow grain size and soot concentration in snow, a mixture of hexagonal columns and plates with rough surfaces.
Overall, while it is clear that spheres are not an ideal choice for modelling the SSPs of snow, it is less clear which non-spherical model should be used. Kokhanovsky and Zege (2004) noted that the final decision of the shape model should be made when in situ phase function measurements for snow become available. The present paper makes a step towards this direction. We employ angular scattering measurements for blowing snow performed with a polar nephelometer  during the CLimate IMpacts of Short-Lived pollutants In the Polar region (CLIMSLIP) campaign at Ny Ålesund, Svalbard (Guyot et al., 2015), to construct a reference phase function for snow grains at the wavelength λ = 0.80 µm. This phase function is used to select a new shape model for snow, an "optimized habit combination" (OHC) consisting of severely rough (SR) droxtals, aggregates of SR plates and strongly distorted Koch fractals. The SSPs for the OHC are then computed as a function of wavelength and snow grain size, and parameterization equations are developed for the single-scattering co-albedo β = 1 − ω, the asymmetry parameter g and the phase function P 11 . Such parameterizations are of substantial practical significance, as they greatly facilitate the use of the OHC in radiative transfer applications. We are not aware of any such previous parameterizations for representing the snow SSPs.
The outline of this paper is as follows. First, in Sect. 2, the models used to compute the SSPs of Koch fractals, Gaussian spheres and spheres are introduced along with the database of Yang et al. (2013) used for several other shapes. In Sect. 3, the reference phase function for snow is constructed. In Sect. 4, several shape models are compared in terms of their ability to reproduce the reference phase function, and the OHC is selected. In Sect. 5, the SSPs for the OHC are computed as a function of wavelength and snow grain size, and in Sect. 6 parameterization equations are developed. In Sect. 7, the snow SSP parameterization is applied to radiative transfer computations, and comparisons are made to spheres and Koch fractals. Finally, a summary is given in Sect. 8.

Shape models and single-scattering data
Here, several shape models are considered as candidates for representing the SSPs of snow. These include (1) secondgeneration Koch fractals, (2) Gaussian random spheres, (3) nine different crystal habits in the Yang et al. (2013) single-scattering database and, for comparison, (4) spheres. The snow grains are assumed to consist of pure ice (i.e. no impurities such as black carbon are included). The ice refractive index of Warren and Brandt (2008) is employed.
The SSPs (extinction cross section, single-scattering albedo, phase function and asymmetry parameter) of Koch fractals are simulated using the geometric optics code of Macke (1993) (see also Macke et al., 1996). Both regular and distorted Koch fractals are considered. A regular secondgeneration Koch fractal has 144 equilateral triangular surface elements. Distortion is simulated using a statistical approach, in which for each refraction-reflection event the normal of the crystal surface is tilted randomly around its original direction (Macke et al., 1996). The zenith (azimuth) tilt angle is chosen randomly with equal distribution between [0, θ max ] ([0, 360 • ]), where θ max is defined using a distortion parameter t = θ max /90 • . Values of t = 0 (regular), t = 0.18 (distorted) and t = 0.50 (strongly distorted) are considered. The geometric optics solution consists of ray tracing and diffraction parts, which are combined as in Macke et al. (1996). For diffraction, the Fraunhofer (far-field) approximation is employed. Either 3 million (in Sect. 4) or 1 million (in Sect. 5) rays per case (i.e. crystal size, wavelength and degree of distortion) are used for the ray tracing part. Up to p = 12 ray-surface interactions per initial ray are considered (see Sect. 3A in Macke, 1993).
The SSPs of Gaussian random spheres are computed with the geometric optics code of Muinonen et al. (1996). Details P. Räisänen et al.: Parameterization of single-scattering properties of snow 1279 of the Gaussian random sphere shape model are discussed e.g. in Nousiainen and McFarquhar (2004). The shape of the particles is described in terms of three parameters: the relative SD of radius σ , the power-law index ν in the Legendre polynomial expansion of the correlation function of radius (the weight of the nth degree Legendre polynomial P n being c n ∝ n −ν ), and the degree of truncation n max for this polynomial expansion. In broad terms, increasing σ increases the large-scale non-sphericity of the particle, while decreasing ν and increasing n max adds small-scale structure to the particle shape. Four values were considered for σ (0.15, 0.20, 0.25 and 0.30), four for ν (1.5, 2.0, 2.5 and 3.0) and three for n max (15, 25 and 35), which yields 48 parameter combinations. A total of 1 million rays with 1000 realizations of particle shape per case were employed in the ray tracing computations. Diffraction was computed by applying the Fraunhofer approximation to equivalent cross-section spheres.
Recently, Yang et al. (2013) published a comprehensive library of SSPs of non-spherical ice crystals, for wavelengths ranging from the ultraviolet to the far infrared, and for particle maximum dimensions d max ranging from 2 to 10 000 µm. The library is based on the Amsterdam discrete dipole approximation (Yurkin et al., 2007) for small particles (size parameter smaller than about 20) and improved geometric optics (Yang and Liou, 1998;Bi et al., 2009) for large particles. Here, single-scattering properties for nine ice particle habits in the Yang et al. (2013) database are used: droxtals, solid and hollow hexagonal columns, aggregates of 8 columns, plates, aggregates of 5 and 10 plates, and solid and hollow bullet rosettes. For each habit, the SSPs are provided for three degrees of particle surface roughness: completely smooth (CS), moderately rough (MR) and severely rough. The effect of roughness is simulated in a way that closely resembles the treatment of distortion for Koch fractals: the surface slope is distorted randomly for each incident ray, assuming a normal distribution of local slope variations with a SD of 0, 0.03 and 0.50 for the CS, MR and SR particles respectively in Eq. (1) of Yang et al. (2013). In fact, this approach does not represent any specific roughness characteristics but attempts instead to mimic the effects on SSPs due to non-pristine crystal characteristics in general (both roughness effects and irregularities).
For comparison, results are also shown for spheres. The SSPs of spheres are computed using a Lorenz-Mie code (de Rooij and van der Stap, 1984;Mischenko et al., 1999).

Observation-based phase function for blowing snow
We employ as a reference an observation-based phase function for blowing snow. The reference phase function was derived from ground-based measurements conducted during the CLIMSLIP field campaign at Ny Ålesund, Svalbard (Guyot et al., 2015), on 23 and 31 March 2012. The blowing snow case on 23 March was preceded by heavy snow-fall on 22 March, ending during the night of 23 March. The last snowfall before 31 March blowing snow case occurred on 29 March. Consequently, the case of 23 March represents essentially new snow, while on 31 March some snow metamorphism had occurred, and the snowpack was probably denser (although snow density was not measured). The near-surface air temperature ranged from −5 to −9 • C during the 23 March case and from −18 to −20 • C during 31 March. Correspondingly, the wind speeds ranged from 1 to 9 m s −1 on 23 March (median value 4 m s −1 ) and from 5 to 8 m s −1 on 23 March (median value 7 m s −1 ). Mainly cloudy conditions prevailed on 23 March, while 31 March was cloud free. The phase functions discussed below are averages over the entire blowing snow events, which lasted for approximately 10 h (8-18 UTC) on 23 March and 12 h (12-24 UTC) on 31 March.
The angular scattering coefficient (θ s ) [µm −1 sr −1 ] of blowing snow was measured with a polar nephelometer (PN; Gayet et al., 1997;Crépel et al., 1997) on 23 and 31 March 2012 at 31 scattering angles in the 15 • ≤ θ s ≤ 162 • range at a nominal wavelength of λ = 0.80 µm. The corresponding phase function P 11 (θ s ) was obtained by normalizing (θ s ) by the volume extinction coefficient σ ext : Here σ ext was estimated from the PN data following Gayet et al. (2002), with a quoted accuracy of 25 %. The derived phase functions are shown in Fig. 1a. There are only minor differences between the 23 and 31 March cases. In both cases P 11 decreases sharply from 15 to 50 • , then more gradually until 127 • . At larger scattering angles P 11 increases slightly with a local maximum around 145 • (discussed below). Hereafter, the average over the two cases is used as a reference for the modelled phase functions: In Fig. 1b, P ref 11 is compared with three other phase functions: a non-precipitating cirrus case over Southern France in the CIRRUS'98 experiment (Durand et al., 1998) Jourdan et al., 2003) and two phase functions for glaciated parts of nimbostratus over Svalbard in the ASTAR 2004 experiment, corresponding to clusters 6 and 7 in Jourdan et al. (2010). These phase functions were derived from raw PN data using a statistical inversion scheme (Jourdan et al., 2003(Jourdan et al., , 2010. Perhaps as expected, the blowing snow phase function P ref 11 is generally closer to the glaciated nimbostratus phase functions than to the cirrus phase function. In particular, at sideward angles between roughly 55 and 135 • , P ref 11 falls mostly between the two nimbostratus phase functions, while the cirrus phase function exhibits somewhat smaller values. The smallest P 11 in the cirrus and nimbostratus cases occurs at θ s = 120 • , as compared with θ s = 127 • for P ref 11 . All four phase functions then increase until θ s ≈ 140 • , after which the nimbostratus and cirrus phase functions become quite flat. In contrast, P ref 11 shows a local maximum around θ s ≈ 145 • . The origin of the maximum at θ s ≈ 145 • is not clear. While it may, in principle, be caused by scattering by snow grains, this feature is neither captured by any of the particle shapes considered in this study nor present in phase functions measured for laboratory-generated ice crystals in Barkey et al. (2002) and Smith et al. (2015). Rather, it falls between the icebow peak for spherical ice particles near 135 • and a maximum seen for many pristine hexagonal shapes at 150-155 • (see Fig. 3). Curiously, this feature coincides with the scattering maximum of small water droplets with a ∼ 10 µm diameter at 140-145 • . However, water droplets seem like an implausible explanation, since the conditions at the measurement site were subsaturated with respect to liquid water (the relative humidity being roughly 92-95 % on 23 March and 79-87 % on 31 March), and especially the 31 March case was quite cold. Yet the 145 • feature is clearly visible in the measured phase function in both cases. Finally, we cannot discount the possibility that inaccuracy in the PN angular scattering measurements influences this feature. Shcherbakov et al. (2006) report relative accuracy of scattered intensities of 3-5 % between 15 and 141 • , but degrading to 30 % for 162 • , for an experimental setup with low extinction. Thus the phase function derived from the PN measurements is, overall, less reliable near the backscattering direction than in near-forward and side-scattering directions.
Whether the phase function feature at 145 • is an artifact or a real feature caused by scattering by snow should be resolved through further measurements, preferably using some alternative technique. However, in either case, it has only a small impact on the snow SSP parameterizations derived in this paper. This detail cannot be captured by any of the shape models considered, so it is not present in the parameterized phase functions. Its influence on the asymmetry parameter is also modest. Even a complete elimination of the maximum by linear interpolation of P ref 11 between the minima at 127 and 155 • would increase g by only ≈ 0.007.
The size distribution of blowing snow was measured with the Cloud Particle Imager (CPI) instrument (Lawson et al., 2001). The CPI registers particle images on a solid state, one million pixel digital charge-coupled device (CCD) camera by freezing the motion of the particle using a 40 ns pulsed, high power laser diode. Each pixel in the CCD camera array has an equivalent size in the sample area of 2.3 µm. In the present study, the minimum size for the CPI's region of interest is set up to 10 pixels. Therefore particles with sizes ranging approximately from 25 µm to 2 mm are imaged. Figure 2a shows examples of particles imaged by the CPI on 31 March 2012. While some needle-shaped crystals can be spotted, many of the particles are irregular, which also applies to the 23 March 2012 case. It is also noted that many of the particles show rounded edges, possibly related to sublimation during snow metamorphosis. Size distributions derived from the CPI data are shown in Fig. 2b. A lognormal distribution was fitted to the data (averaged over the 23 and 31 March cases): Here, d p is the projected-area equivalent diameter of the particles, d p,0 = 187 µm is the median diameter and σ g = 1.62 the geometric SD. This size distribution was used for all shape models when comparing the modelled phase functions with P ref 11 . Since absorption is weak at λ = 0.80 µm and the particles are much larger than the wavelength, the modelled P 11 is only weakly sensitive to the size distribution employed if the shape of the snow grains is independent of size. This holds true for spheres, Gaussian spheres, Koch fractals, droxtals and the three aggregate habits in the Yang et al. (2013) database. However, for solid and hollow hexagonal columns, plates as well as solid and hollow bullet rosettes, the crystal geometry is a function of size, with some influence on P 11 (see end of Sect. 4 for more discussion).

Selecting a shape model for snow optics
The purpose of this section is to select a shape model of snow for use in Sects. 5 to 7. The phase function for blowing snow from the CLIMSLIP campaign, as defined by Eq.
(2), is used as a reference. It is emphasized that the approach is deliberately pragmatic: we do not attempt to model the scattering based on the shapes of the observed snow grains; rather we try to develop an equivalent microphysical model for representing the SSPs. Previously, the choice of Koch fractals for approximating the scattering by snow (Kokhanovsky and Zege, 2004) was likewise based on phase function data only. Furthermore, our approach is conceptually analogous to the widely used practice of modelling the SSPs of irregular dust particles. Instead of considering the actual dust particle shapes, shape distributions of spheroids are used operationally in a variety of applications (Dubovik et al., 2006(Dubovik et al., , 2011Levy et al., 2007), as they have been found to reasonably mimic scattering by dust. In contrast, current state-of-the-art models for ice cloud SSPs include ice crystal habit distributions parameterized as a function of crystal size, based on in situ microphysical observations (Baum et al., 2005(Baum et al., , 2011Hong et al., 2009). In principle, it would be desirable to use this approach also for snow to provide a more direct link between the actual snow grain shapes and those assumed in the parameterization and to account for changes in snow grain shape with size, which we currently neglect. This would require, first, the analysis and subsequent parameterization of snow grain shape distributions as a function of size and, second, the computation and parameterization of the respective SSPs. The main reason why we have not attempted this approach in the current work is that a very large fraction of the particles in blowing snow (and snow on ground) are irregular, more than 80 % according to manual classification of CPI images (Guyot et al., 2015) (see also (4) To start with, the phase function for single crystal shapes is compared with P ref 11 in Fig. 3. To be consistent with the CLIMSLIP observations, the phase function is computed at λ = 0.80 µm, and it is integrated over the size distribution defined by Eq. (3). Several points can be noted.
First, unsurprisingly, the phase function for spheres agrees poorly with the observations (Fig. 3a). In particular, sideward scattering is underestimated drastically, and there is a strong icebow peak at θ s = 134 • , which is not seen in P ref 11 . Second, for second-generation Koch fractals (Fig. 3b), the agreement with P ref 11 is considerably better than for spheres. The main features of the phase function are similar for regular and distorted Koch fractals. However, the regular Koch fractal's phase function exhibits several sharp features specific to the tetrahedral geometry which are not observed in P ref 11 . The distorted Koch fractals' versions are more consistent with the measurements even though marked deviations from P ref 11 are still present. Scattering is underestimated between 15 and 30 • and overestimated between 45 and 100 • .  Yang et al. (2013) database. For each habit, the phase function was averaged over the size distribution defined by Eq. (3). In the figure legends, the two numbers in parentheses give the asymmetry parameter and the cost function defined by Eq. (4) respectively. For the Gaussian spheres in (c), the notation indicates the shape parameters (e.g. for 0.15_3.0, σ = 0.15 and ν = 3.0; n max was fixed at 15). For the Yang et al. (2013) habits in (d-l), CS, MR and SR refer to particles with completely smooth surface, moderate surface roughness and severe surface roughness respectively. Also, the gradient of P 11 in the backscattering hemisphere is consistently negative, while P ref 11 rather increases slightly between 127 and 162 • . Overestimated sideward scattering by Koch fractals has been previously noted in the context of cirrus clouds (Francis et al., 1999) and in a comparison with a measured phase function for laboratory-generated ice crystals ( Fig. 3 in Kokhanovsky and Zege, 2004).
Third, for Gaussian spheres, the level of agreement with P ref 11 depends on the shape parameters chosen. Four cases out of the 48 considered are shown in Fig. 3c (for all of these n max = 15, but the general features for n max = 25 and n max = 35 are similar). For example, for the parameter values σ = 0.15 and ν = 3.0, which are close to those estimated from shape analysis of small quasi-spherical ice crystals in cirrus clouds in Nousiainen et al. (2011), the deviations from P ref 11 are substantial. The phase function features undesirable large-scale oscillations and, in particular, scattering at θ s ≈ 45-75 • is underestimated substantially. Best agreement with P ref 11 is obtained in the case σ = 0.30, ν = 0.15, which features both pronounced large-scale non-sphericity and small-scale structure in the particle shape. The sideward scattering is overestimated (mainly between 70 and 100 • ), but the cost function (0.163) is clearly smaller than that for distorted Koch fractals (0.284) and is, in fact, the smallest among all single-habit shape models considered.
Fourth, regarding the habits in the Yang et al. (2013) database ( Fig. 3d-l), both visual inspection and the cost function values indicate that the agreement with P ref 11 improves with increasing particle surface roughness. While completely smooth and, in many cases, moderately rough particles exhibit halo peaks, for severely rough particles the phase function is quite smooth and featureless, as is P ref 11 . It is further seen that, in general, increasing the roughness increases sideward scattering and reduces the asymmetry parameter. While none of the habits considered provides perfect agreement with P ref 11 , the cost function is smallest for the aggregate of eight columns (0.172).
Since none of the individual shape models agrees fully satisfactorily with P ref 11 , we considered combinations of two or three shapes. We thus use where n = 2 or n = 3 is the number of shapes in a combination and P j 11 is the phase function for shape j , integrated over the size distribution (Eq. 3) for each shape separately. Thus, the potential dependence of snow grain shapes on their size is not considered here. For each combination of shapes considered, the optimal weight factors w j were searched by minimizing the cost function (Eq. 4), subject to the conditions that all w j are non-negative and their sum equals 1. Since pristine particles and even moderately rough particles feature halo peaks (or an icebow peak in the case of spheres), which are absent in P ref 11 , the following groups of habits are considered: distorted Koch fractals, Gaussian spheres and severely rough particles in the Yang et al. (2013) database. Figure 4 illustrates a comparison with P ref 11 for three singlehabit cases ( Fig. 4a and d) (the best Koch fractal case, the best Gaussian sphere case and the best case with Yang et al. Table 1. List of best three-habit combinations. w 1 , w 2 and w 3 are the weights (i.e. fractional contributions to projected area) of each habit, "cost" is the cost function, g the asymmetry parameter and ξ a non-dimensional absorption parameter defined by Eq. (8). SR refers to severely rough particles, and t is the distortion parameter for second-generation Koch fractals. The "optimized habit combination" (OHC) is highlighted with italic font. (2013) particles), the best three two-habit cases ( Fig. 4b and e) and the best three three-habit cases ( Fig. 4c and f) as defined in terms of the cost function. As expected, the agreement of P model 11 with P ref 11 improves with increasing number of shapes in the combination. The best three-habit cases follow P ref 11 quite faithfully, though slightly underestimating P ref 11 in near-forward directions and not capturing the details of P ref 11 near θ s = 145 • . Furthermore, it is seen that the best three-habit combinations produce nearly identical P 11 , agreeing even better with each other than with P ref 11 . The relationship between the asymmetry parameter g and the cost function is considered in Fig. 5, where all singlehabit cases and combinations of two or three habits are included. While high values of cost function can occur at any g, the lowest values (< 0.10) always occur for three-habit combinations with 0.775 < g < 0.78. This supports a best estimate of g ≈ 0.78 for snow at λ = 0.80 µm, of course subject to the assumption that the measurements for blowing snow used to construct P ref 11 are also representative of snow on the ground.
The three-habit combinations with cost function below 0.1 are listed in Table 1. All of them include SR droxtals and either strongly distorted (t = 0.50) or distorted (t = 0.18) Koch fractals, but the third habit included in the combinations varies from case to case. The differences in cost function and asymmetry parameter between the best habit combinations are very small, which makes the choice of a single "best" habit combination for representing the SSPs of snow somewhat arbitrary. For further use in representing the SSPs as a function of wavelength and size, we select the following habit combination: 36 % of SR droxtals, 26 % of aggregates of 10 SR plates and 38 % strongly distorted secondgeneration Koch fractals (t = 0.50), where the weights refer to fractional contributions to the projected area. This habit combination is represented with a blue line in Fig. 4c and f and is marked with an arrow in Fig. 5. Hereafter, this habit combination will be referred to as the "optimized habit com- Figure 5. A scatter plot of asymmetry parameter vs. cost function (Eq. 4) for single habits (black dots), combinations of two habits (red dots) and combinations of three habits (blue dots). The "optimized habit combination" selected for parameterization of snow single-scattering properties is marked with an arrow. Note that some single-habit cases fall outside the range plotted here. These include spheres for which cost = 1.90 and g = 0.892. bination". The primary reason why we selected this OHC rather than either of the first two habit combinations in Table 1, which have a marginally lower cost function, is that these habit combinations include either hollow columns or The Cryosphere, 9, 1277-1301, 2015 www.the-cryosphere.net/9/1277/2015/ bullet rosettes. For these habits (unlike aggregates of plates), the particle geometry assumed in the Yang et al. (2013) database depends on particle size, with the aspect ratio of the crystals increasing with their length. However, due to snow metamorphosis on ground, size-shape relationships based on crystal growth in ice clouds are most likely not valid for snow. Therefore, we considered it better to use a crystal geometry that is independent of size. This also helps to keep the SSP parameterization simpler.

Snow single-scattering properties as a function of size and wavelength
The SSPs, including the extinction efficiency Q ext , singlescattering co-albedo β, asymmetry parameter g and scattering phase function P 11 (θ s ) were determined for the OHC for 140 wavelengths between 0.199 and 3 µm and for 48 particle sizes between 10 and 2000 µm. Here, the size is defined as the volume-to-projected area equivalent radius r vp = 0.75 V /P . As stated above, the OHC consists of SR droxtals, aggregates of 10 SR plates and strongly distorted Koch fractals. The SSPs for droxtals and aggregates of plates were taken from the Yang et al. (2013) database (interpolated to fixed values of r vp ) while those of Koch fractals were computed using the geometric optics code of Macke (1993), as explained in Sect. 2. Four caveats should be noted: 1. Due to problems associated with the truncation of numerical results to a finite number of digits (P. Yang, personal communication, 2013), the values of β in the Yang et al. (2013) database are unreliable in cases of very weak absorption. To circumvent this issue, it was assumed that in cases of weak absorption (β < 0.001 for Koch fractals), the values for droxtals and aggregates of plates may be approximated as β aggregate λ, r vp = 0.932β fractal λ, r vp .
Here the scaling factors were determined as β droxtal /β fractal and β aggregate /β fractal , where the overbar refers to averages over the cases in which 0.001 < β fractal < 0.01 and the size parameter x = 2π r vp /λ > 100. 3. The SSPs for Koch fractals were computed using a geometric optics code, which means that the accuracy deteriorates somewhat in cases with smaller size parameters (typically for x < 100). This issue pertains mainly to small snow grains at near-IR wavelengths (e.g. for λ = 2.5 µm, x = 100 corresponds to r vp ≈ 40 µm).
4. Lastly but importantly, since the OHC was selected based on measurements at a single wavelength λ = 0.80 µm for only two cases, there is no guarantee that it represents the snow SSPs equally well at other wavelengths, or for all snow grain sizes. Figure 6 compares wavelength-dependent SSPs for the OHC with those for two shape assumptions previously used in modelling snow optics: spheres and Koch fractals (distorted Koch fractals with t = 0.18 were selected for this comparison; this is close though not identical to the shape assumption used by Kokhanovsky et al., 2011). Two monodisperse cases are considered, with r vp = 50 µm and r vp = 1000 µm. For all three habits, the asymmetry parameter g (Fig. 6a) and the single-scattering co-albedo β (Fig. 6b) show well-known dependencies on particle size and wavelength. Thus, g is largely independent of both λ and r vp in the visible region where β is very small. In the near-IR region, β increases with increasing imaginary part m i of the refractive index and with increasing particle size. With increasing β, the fractional contribution of diffraction to the phase function increases, which results in larger values of g (e.g. Macke et al., 1996). The most striking differences between the three shape assumptions occur for the asymmetry parameter, especially in the visible region, where g ≈ 0.89 for spheres, g ≈ 0.74 for distorted Koch fractals and g ≈ 0.77-0.78 for the OHC. The values of β for the OHC are also intermediate between the two single-shape cases: larger than those for spheres (except for r vp = 1000 µm at the strongly absorbing wavelengths λ > 1.4 µm) but slightly smaller than those for distorted Koch fractals. The implications of these differences for snow albedo are considered in Sect. 7.
While the co-albedo values in Fig. 6b are strongly wavelength dependent through m i , the effects of shape on absorption can be distinguished more clearly by considering the non-dimensional absorption parameter (Kokhanovsky and Zege, 2004;Kokhanovsky, 2013) where C abs is the absorption cross section, Q ext the extinction efficiency, P the projected area and V the particle volume, and γ = 4π m i /λ, where m i is the imaginary part of ice refractive index. Figure 6c displays ξ at the wavelengths λ = 0.199-1.4 µm, where absorption by snow is relatively weak. Consistent with the co-albedo values (Fig. 6b) and previous studies (e.g. Kokhanovsky and Nauss, 2005), Fig. 6c indicates that absorption is generally stronger for nonwww.the-cryosphere.net/9/1277/2015/ The Cryosphere, 9, 1277-1301, 2015 spherical than spherical particles for the same r vp . The difference is particularly clear in the visible region, where ξ ≤ 1.3 for spheres (except for some spikes that occur in the Mie solution especially for r vp = 50 µm), ≈ 1.7 for the Koch fractals and slightly over 1.6 for the OHC. At wavelengths beyond λ = 1.0 µm, ξ tends to decrease especially for the larger particle size r vp = 1000 µm considered, as absorption no longer increases linearly with m i . Furthermore, in the UV region, Koch fractals and the OHC show a distinct increase in ξ with decreasing wavelength. This is related to the corresponding increase of the real part of the refractive index m r . Interestingly, it is found that for these shape assumptions, absorption scales linearly with m 2 r ; furthermore, for Koch fractals ξ/m 2 r ≈ 1 when absorption is weak (Fig. 6d). For spheres, the dependence of ξ on m r is weaker. Equation (4) in Bohren and Nevitt (1983) provides the absorption efficiency of weakly absorbing spheres in the limit of geometric optics, which can be rewritten in terms of ξ as For r vp = 1000 µm, ξ for spheres follows this approximation closely until λ ≈ 1.0 µm ( Fig. 6c and d). However, it appears that for Koch fractals, only the first term should be included. It should be noted that ξ for the OHC is not independent of that for Koch fractals (due to the scaling of co-albedo in Eqs. (6) and (7)). However, we found that ξ also scales linearly with m 2 r for Gaussian spheres (this was tested for σ = 0.17, ν = 2.9, n max = 15), suggesting that this might apply more generally to complex non-spherical particles.
Finally, it should be recalled that our choice of the OHC was based on phase function observations at the wavelength λ = 0.80 µm. At this wavelength, absorption is so weak that it has very little impact on the phase function. Therefore, these observations cannot be used to constrain absorption by snow. In spite of this, we think it is worth providing a co-albedo parameterization based on the OHC (Eq. 11 in Sect. 6.2). The reason for this is that snow grains are distinctly non-spherical, and for non-spherical particles, ξ and β are, in general, systematically larger than those for spheres, as demonstrated by Fig. 6. In fact, considering the wavelength λ = 0.80 µm, the values of ξ integrated over the size distribution defined by Eq. (3) are, for the large majority of the non-spherical shapes considered, between 1.55 and 1.75, the value for the OHC being ξ = 1.62 (Table 1). The corresponding value for spheres is substantially lower: ξ = 1 .29. Thus, while we cannot constrain ξ or β precisely, it is very likely that the actual values for snow exceed those for spheres.

Parameterizations for the single-scattering properties of snow
In this section, parameterization equations are provided for the computation of snow SSPs (extinction efficiency Q ext , single-scattering co-albedo β, asymmetry parameter g and scattering phase function P 11 (θ s )) for the OHC discussed above. The parameterizations are provided for the size range r vp = 10-2000 µm and wavelength range λ = 0.199-2.70 µm. They are expressed in terms of the size parameter x and real and imaginary parts of refractive index (m r and m i ). Here, the size parameter defined with respect to the volume-toprojected area equivalent radius is used: For the OHC, the size parameter defined with respect to the projected area is x p ≈ 1.535 x vp .

Extinction efficiency
The extinction efficiency Q ext for the OHC is displayed in Fig. 7. For most of the wavelength and size region considered, Q ext is within 1 % of the asymptotic value Q ext = 2 for particles that are large compared to the wavelength. Note that the deviations from Q ext = 2 are probably somewhat underestimated because the OHC includes Koch fractals, for which Q ext ≡ 2 due to the use of geometric optics. For simplicity, we assume this value in our parameterization, while acknowledging that the actual value tends to be slightly higher especially for small snow grains in the near-IR region.

Single-scattering co-albedo
The single-scattering co-albedo is parameterized as where the size parameter for absorption is defined as Figure 7. Extinction efficiency Q ext for the optimized habit combination as a function of wavelength (λ) and volume-to-projected area equivalent radius (r vp ).
The general form of this parameterization was inspired by the ice crystal optics parameterization of van Diedenhoven et al. (2014); however, our definition of x abs differs from theirs in that the factor m 2 r is included based on the findings of Fig. 6c and d. The performance of this parameterization is evaluated in Fig. 8a and c. In Fig. 8a, the parameterized values (shown with contours) follow extremely well the reference values computed for the OHC (shading). The relative errors β/β are mostly below 1 %; errors larger than 3 % (and locally even > 10 %) occur only for small snow grains (r vp < 50 µm) at wavelengths λ > 1.2 µm. The rms value of the relative errors (computed over 125 values of λ ∈ [0.199, 2.7 µm] and 48 roughly logarithmically spaced values of r vp ∈ [10, 2000 µm]) is 1.4 %.

Asymmetry parameter
The asymmetry parameter is parameterized as where the parameter values were determined by trial and error, with the aim of minimizing the rms error in g. The form of this parameterization reflects how g decreases with increasing m r , increases with increasing absorption (i.e. increasing co-albedo β) and increases slightly with increasing size parameter x vp even at non-absorbing wavelengths, in part because the diffraction peak becomes narrower. In practice, the co-albedo β plays the most important role (cf. van Diedenhoven et al., 2014), which explains the general increase of g with increasing r vp in the near-IR region (Fig. 8b). The parameterized values of g (shown with contours in Fig. 8b) follow the reference values (shading) very well. Note that when producing these results, parameterized rather than exact β was used in Eq. (13). The differences from the reference are mostly below 0.001 at the weakly ab- sorbing wavelengths up to λ = 1.4 µm, and while larger differences up to |g| = 0.007 occur at the strongly absorbing wavelengths (Fig. 8d), the overall rms error is only 0.0019.

Phase function
The phase function parameterization consists of three terms, which represent contributions due to diffraction, due to the ray tracing part and a residual that corrects for errors made in approximating the former two parts. The weight factors for diffraction w diff and ray tracing w ray are given by where the latter form assumes Q ext = 2 (e.g. Macke et al., 1996). It should be noted that in practice, the division of the phase function expressed by Eq. (14) is conceptual rather than rigorous. The fitting was based on the total phase function rather than the diffraction and ray tracing parts separately, as these two parts are not separated in the Yang et al. (2013) database. The general aim of the fitting was to minimize the rms errors in ln P 11 .
For diffraction, the HG phase function (Henyey and Greenstein, 1941) is used: The HG phase function is given by and the asymmetry parameter g diff is approximated as a rough approximation and clearly not ideal for studies of very near-forward scattering, but it serves well the current purpose. On one hand, it improves the accuracy compared to the assumption of a delta spike, and on the other hand, the HG phase function has a very simple Legendre expansion P HG (g, θ s ) = ∞ n=0 (2n + 1)g n P n (cos θ s ) , where P n denotes the nth order Legendre polynomial. This facilitates greatly the use of P HG in radiative transfer models such as DISORT (Stamnes et al., 1988). The phase function for the ray tracing part is approximated as where the latter term 1 − w 1 is intended to emulate the nearly flat behaviour of P 11 in the near-backward scattering directions. The weight factor for the HG part is parameterized as where g ray is the asymmetry parameter for the ray tracing (i.e. non-diffraction) part. It is derived from the condition g = w diff g diff + w ray g ray , which yields The total asymmetry parameter g is computed using Eq. (13) above. Finally, the asymmetry parameter g 1 needed in Eq. (21) is While the sum of the first two terms of Eq. (14) already provides a reasonably good approximation of the phase function (see below), the fit can be further improved by introducing the residual P resid , which is represented as a Legendre series. It turns out that, except for cases with strong absorption, a series including terms only up to n = 6 yields very good results P resid (θ s ) = 6 n=0 (2n + 1)a n P n (cos θ s ) , provided that δ-M scaling (Wiscombe, 1977) is applied, with a truncated fraction f = a 6 . Thus, (2n + 1) (a n − a 6 ) P n (cos θ s ) , (26) where δ is Dirac's delta function. What remains to be parameterized, then, are the coefficients a 0 . . . a 6 . A rough but useful approximation is to express them as a simple function of the co-albedo β and the asymmetry parameter g: a n = c 1n + c 2n β + c 3n g + c 4n βg.
The parameterization coefficients c mn were determined by minimizing the rms errors of a n with the LAPACK subroutine DGELS, and they are given in Table 2. Note specifically that the coefficients c m0 and c m1 are all 0. The formulation of P diff and P ray ensures that the phase function (Eq. 14) is correctly normalized and that its asymmetry parameter is consistent with Eq. (13) even without considering P resid ; therefore a 0 = a 1 = 0. Equivalently, the Legendre expansion may be replaced by an ordinary polynomial. This yields Here, the coefficients d mn were obtained directly based on the coefficients c mn in Eq. (27) by writing out the Legendre polynomials in Eq. (26). Their numerical values are given in Table 3. In summary, the phase function parameterization reads P 11 (θ s ) = w diff P HG (g diff , θ s ) + w ray w 1 P HG (g 1 , θ s ) where P resid (θ s ) is given by Eq. (26) or, equivalently, by Eq. (28). Finally, it is worth noting how this parameterization can be used in DISORT, when applying a "δ-NSTR-stream" approximation for radiative transfer, NSTR being the number of streams. In this case, DISORT assumes by default a truncation factor f = a NSTR . If NSTR > 6, the Legendre expansion for P resid in Eq. (26) should be formally extended to n = NSTR, with a n = a 6 for n = 7 . . . NSTR. Thus the Legendre coefficients input to DISORT become for n = 0 w diff g n diff + w ray w 1 g n 1 + a n , for 1 ≤ n ≤ 6 where we have utilized the Legendre expansion of the HG phase function in Eq. (20).
To provide a compact view of how the phase function parameterization performs, we define, analogously to Eq.   is the parameterized phase function and P OHC 11 is the reference value, defined here as the "exact" phase function computed for the OHC. Figure 9a shows the cost function for the full phase function parameterization, and Fig. 9b shows that for a simpler parameterization that includes only the first two terms of Eq. (14) (i.e. P resid is excluded). Note that the parameterized phase function is computed here using parameterized (rather than exact) values of Q ext , β and g.
Most importantly, Fig. 9a shows that in a large part of the wavelength and size domain, the accuracy of the full parameterization is very high, with cost function values ≤ 0.03. This corresponds to a typical relative accuracy of 3 % in the computed phase function, as compared with the reference values for the OHC. The primary exception is that substantially larger errors occur for large snow grains at the strongly absorbing wavelengths in the near-IR region. In broad terms, the accuracy starts to degrade appreciably when β > 0.3, that is, in cases in which snow reflectance is quite low (β = 0.3 corresponds roughly to a spherical albedo of 0.03 for an optically thick snow layer). At the largest wavelengths considered (λ > 2.5 µm), somewhat larger values of the cost function also occur for smaller values of r vp and β. The cost function for the simplified parameterization (Fig. 9b) shows mainly the same qualitative features as the full parameter- ization in Fig. 9a; however, the cost function values in the weakly absorbing cases are ≈ 0.07, in contrast to the values of ≈ 0.03 for the full parameterization. Figure 10 displays examples of phase function for nine combinations of λ and r vp . In the weakly absorbing cases in Fig. 10a-c, and also at the more strongly absorbing wavelength λ = 1.50 µm for r vp = 10 µm and r vp = 100 µm ( Fig. 10d and e), the full parameterization follows extremely well the reference phase function computed for the OHC, to the extent that the curves are almost indistinguishable from each other. Even at λ = 2.00 µm, the deviations from the reference are generally small in the cases with relatively small snow grains (r vp = 10 µm and r vp = 100 µm; Fig. 10g and h), although backward scattering is slightly overestimated in the latter case. In contrast, in cases with very strong absorption and large snow grains (r vp = 1000 µm for λ = 1.50 µm and λ = 2.00 µm in Fig. 10f and i) there are more substantial deviations from the reference. Here, the parameterized phase function is generally underestimated in the backscattering hemisphere and overestimated at θ s < 30 • especially for λ = 2.00 µm, r vp = 1000 µm. Furthermore, the Legendre expansion in P resid leads to oscillations in the backscattering hemisphere which do not occur in the reference phase function. Again, it should be noted that the largest errors occur Figure 10. Examples of the reference phase function computed for the OHC (black lines) and of the parameterized phase function for the full parameterization (red lines), the simplified parameterization without the term P resid in Eq. (14) (blue lines) and the Henyey-Greenstein phase function with asymmetry parameter defined by Eq. (13) (dashed green lines) for nine combinations of wavelength λ and volume-to-projected area equivalent radius r vp . For reference, the values of single-scattering co-albedo β, asymmetry parameter g and cost functions for the full parameterization (cost1), for the simplified parameterization (cost2) and for the Henyey-Greenstein phase function (cost3) are listed in each panel.
in cases in which snow is very "dark": the spherical albedo corresponding to the cases in Fig. 10f and i is only ∼ 0.005.
In many respects, the simplified parameterization (i.e. without P resid ) produces quite similar phase functions as the full parameterization. Two differences can be noted. First, the simplified parameterization does not capture the slight increase in phase function at angles larger than θ s ≈ 120-130 • , which is present in the reference and full parameterization phase functions and was also suggested by the CLIMSLIP data for blowing snow at λ = 0.80 µm, along with the other phase functions in Fig. 1b. Second, in the cases with very strong absorption ( Fig. 10f and i) the simplified phase function avoids the oscillations seen in the full parameterization.
The utility of providing a phase function parameterization is further demonstrated by showing in Fig. 10, for comparison, the HG phase function computed using the asymmetry parameter from Eq. (13). The differences from the reference phase function are systematic. The scattering in the diffraction peak is underestimated (although this is not properly seen from Fig. 10), but otherwise forward scattering is overestimated until a scattering angle of ≈ 35-80 • , depending on the case. Conversely, at sideward and backscattering angles, scattering is underestimated. Consequently, the cost function values for the HG phase function given in Fig. 10  substantially exceed those for both the full and simplified phase function parameterizations.

Radiative transfer applications
In this section, we consider the impact of snow optics assumptions on snow spectral albedo A and reflected radiances L ↑ .The purpose is, on one hand, to evaluate the accuracy of the proposed snow SSP parameterization and, on the other hand, to compare the results obtained with three shape assumptions: spheres, second-generation Koch fractals (distorted with t = 0.18) and the OHC proposed here. Throughout this section, the results for the OHC are used as the reference, although it is clear that they cannot be considered an absolute benchmark for scattering by snow. The radiative transfer computations were performed with DIS-ORT (with 32 streams, δ-M scaling included), assuming an optically thick (i.e. semi-infinite) layer of pure snow with a monodisperse size distribution.
Like most other solar radiative transfer studies involving snow, close-packed effects are ignored in the calculations. It has been shown by Kokhanovsky (1998) that, at least as a first approximation, they do not have a pronounced impact on the snow reflectance.
First, snow albedo as a function of λ and r vp is considered in Fig. 11. Direct incident radiation with a cosine of zenith angle µ 0 = cos θ 0 = 0.5 is assumed. Figure 11a demonstrates the well-known features of snow albedo: the values are very high in the UV and visible region and decrease with increasing particle size in the near-IR. The results computed using the parameterized snow optical properties Q ext , β, g and P 11 are almost indistinguishable from those obtained using the "exact" optical properties for the OHC. The differences between these two are mostly within 0.002 (Fig. 11b), although larger differences up to 0.02 occur for very small snow grains (r vp ≈ 10-20 µm) at wavelengths with strong absorption by snow (λ > 1.4 µm). These results are only weakly sensitive to the assumed direction of incident radiation. Furthermore, while the parameterized albedo values were computed using the full phase function parameterization, the values for the simplified parameterization (without P resid in Eq. (14)) differed very little from them, mostly by less than 0.001.
For distorted Koch fractals, the albedo values are higher than those for the OHC, but the difference is rather small, at most 0.017 (Fig. 11c). Conversely, for spheres the albedo values are lower, with largest negative differences of −0.08 from the reference (Fig. 11d). This stems from the higher asymmetry parameter of spheres, which is only partly compensated by their lower co-albedo (Fig. 6). To put it another way, for a given albedo A in the near-IR region, a smaller (slightly larger) particle size is required for spheres (for distorted Koch fractals) than for the OHC.
To compare the simulated radiance distributions to the reference, we next consider the root-mean-square error in the logarithm of reflected radiances integrated over the hemisphere: where θ and φ denote the zenith angle and azimuth angle respectively and L ↑ OHC is the radiance in the reference computations for the OHC. Figure 12a-c show LO-GRMSE as a function of particle size and wavelength for the full parameterization for three directions of incident radiation (µ 0 = 0.8, µ 0 = 0.4 and µ 0 = 0.1, corresponding to θ 0 = 36.9 • , θ 0 = 66.4 • and θ 0 = 84.3 • ). For weakly absorbing wavelengths up to λ = 1.4 µm, the performance of the parameterization is extremely good for all particle sizes, with values of LOGRMSE < 0.01 for µ 0 = 0.8 and µ 0 = 0.4 and between 0.01 and 0.02 for µ 0 = 0.1. LOGRMSE ∼ 0.01 implies a typical relative accuracy of ∼ 1 % in the reflected radiances. The accuracy in radiances at weakly absorbing wavelengths is even higher than that in the phase function ( Fig. 9a) because strong multiple scattering diminishes the effect of phase function errors. At wavelengths λ > 1.4 µm, LOGRMSE increases, not only due to larger phase function errors but also because multiple scattering is reduced due to stronger absorption. Even here, LOGRMSE stays mainly below 0.05 for relatively small snow grains (r vp < 100 µm), but substantially larger errors occur in the cases with large and strongly absorbing grains, consistent with the modest accuracy of the phase function parameterization in these cases (Fig. 9a). These errors depend only weakly on µ 0 . It should be noted that the largest relative errors occur in cases where the reflected radiances and radiance errors are small in an absolute sense and probably matter little for practical applications.
Values of LOGRMSE obtained using the simplified phase function parameterization are shown in Fig. 12d-f. Consistent with the phase function errors (cf. Fig. 9a vs. b), the simplified parameterization is slightly less accurate in simulating reflected radiances than the full parameterization except for the most strongly absorbing cases. Nevertheless, the accuracy is quite high for the weakly absorbing cases; LO-GRMSE ranging from ∼ 0.01 (or even less) for µ 0 = 0.8 to ∼ 0.03 for µ 0 = 0.1.
For comparison, Fig. 12g and h show LOGRMSE computed for distorted Koch fractals and spheres (for µ 0 = 0.4 only). Unsurprisingly, LOGRMSE is generally smaller for Koch fractals than for spheres (e.g. 0.05-0.10 in weakly absorbing cases compared to ∼ 0.20 for spheres). In both cases, again excepting large particles at strongly absorbing wavelengths, the values of LOGRMSE are substantially larger than those associated with the snow SSP parameterization. This indicates that in general, numerical fitting errors in the parameterization are a minor issue in comparison with the radiance differences associated with different shape assumptions.
Examples of the angular distribution of reflected radiances are given in Figs. 13 and 14. Here, only a single particle size r vp = 200 µm is considered, and the azimuth angle for incident radiation is φ 0 = 0 • . In Fig. 13, results are shown for three zenith angles of incident radiation, corresponding to µ 0 = 0.8, µ 0 = 0.4 and µ 0 = 0.1, for a single wavelength λ = 0.80 µm. In Fig. 14, three wavelengths are considered (λ = 0.30, 1.40 and 2.20 µm) but for µ 0 = 0.4 only. In each figure, panels a-c display the distribution of reflected radiances in the reference calculations for the OHC, while the remaining panels show the relative differences from the reference for distorted Koch fractals with t = 0.18 (panels d-f), for spheres (g-i), for the Henyey-Greenstein phase function (j-l), for the full snow SSP parameterization (m-o) and for the simpler parameterization without P resid in Eq. (14) (p-r). For brevity, only some main points are discussed.
First, it is seen, consistent with Fig. 12, that in general the radiance distribution for spheres differs more from the reference than the distribution for Koch fractals does. For example, for λ = 0.80 µm and µ 0 = 0.4 both positive and negative differences larger than 50 % occur for spheres (Fig. 13h), while for Koch fractals the differences exceed 10 % only locally (Fig. 13e). Furthermore, in the same case, the radiance errors are < 1 % almost throughout the (θ , φ) domain for the full parameterization (Fig. 13n) and mostly < 2 % even for the simplified parameterization (Fig. 13q). In contrast, when the HG phase function is employed in the calculations, the differences from the reference reach locally 30 and −40 % (Fig. 13k).
Second, while the results noted above for λ = 0.80 µm and µ 0 = 0.4 are also mostly valid for µ 0 = 0.8 and µ 0 = 0.1 and for λ = 0.30, 1.40 and 2.20 µm, some quantitative differences can be noted. When µ 0 decreases from 0.8 to 0.1, the pattern of reflected radiances becomes increasingly non-uniform and more sensitive to both the assumed particle shape and the errors in phase function parameterization. This occurs because the relative role of first-order scattering increases (e.g. Mischenko et al., 1999). For the same reason, the sensitivity of the radiance pattern to the phase function increases with increasing absorption. Thus, while the qualitative features are mostly similar at all wavelengths considered here, the relative differences are generally larger at λ = 1.40 µm and λ = 2.20 µm than at λ = 0.30 µm and λ = 0.80 µm. Espe- cially at the wavelength λ = 2.20 µm, at which snow absorption is quite strong and the albedo for the OHC is only 0.11, the radiance pattern is dominated by first-order scattering and is thus very sensitive to the details of the phase function. In a relative (though not absolute) sense, the errors in parameterized radiances are also somewhat larger than at the other wavelengths considered (Fig. 14o and r).
Third, even at weakly absorbing wavelengths, the role of first-order scattering is clearly discernible: many differences in the pattern of reflected radiances can be traced directly to phase function differences. For example, considering the results for λ = 0.80 µm for both µ 0 = 0.4 and µ 0 = 0.1, we note the following.
-Three regions appear in the radiance differences between distorted Koch fractals and the OHC in Fig. 13e and f. Going from left to right, negative radiance differences occur at large values of θ and small values of φ (roughly for θ > 65 • and φ < 20 • ), followed by a region of positive differences and another region of negative differences (roughly for θ > 40 • , φ > 140 • ). These regions occur because the phase function for Koch fractals is larger than that for the OHC at intermediate scattering angles (29 • ≤ θ s ≤ 134 • ) but smaller in the near-forward and near-backward directions.
-For spheres in Fig. 13h and i, the reflected radiances greatly exceed those for the OHC for roughly θ > 60 • , φ < 40 • because the phase function for spheres is generally larger than that for the OHC for θ s < 54 • . Conversely, at larger θ s the phase function for spheres is (mostly) considerably smaller than that for the OHC. This results in generally smaller reflected radiances for www.the-cryosphere.net/9/1277/2015/ The Cryosphere, 9, 1277-1301, 2015 spheres in most of the (θ , φ) domain with φ > 50 • . As an exception, the icebow feature for spheres at θ s ≈ 135 • results in an arc with larger radiances for spheres than for the OHC.
-For the HG phase function, the pattern of overestimated radiances up to φ ∼ 60 • and underestimated radiances at larger azimuth angles ( Fig. 13k and l) arises because the HG phase function exceeds that for the OHC for θ s < 80 • and falls below it at larger scattering angles (see also Fig. 10).

Summary
In this work, measurements of angular distribution of scattering by blowing snow made during the CLIMSLIP campaign in Svalbard were used to select a shape model for represent-ing the single-scattering properties of snow. An optimized habit combination consisting of SR droxtals, aggregates of SR plates and strongly distorted Koch fractals was selected. The SSPs (extinction efficiency Q ext , single-scattering coalbedo β, asymmetry parameter g and phase function P 11 ) were then computed for the OHC as a function of wavelength and snow grain size. Furthermore, parameterization equations were developed for the SSPs for the wavelength range λ = 0.199-2.7 µm and for snow grain volume-to-projected area equivalent radii r vp = 10-2000 µm. The parameterizations are expressed in terms of the size parameter and real and imaginary parts of refractive index. The relative accuracy of the parameterization, as compared with the reference calculations for the OHC, is very high for the single-scattering coalbedo and the asymmetry parameter. This is also true for the phase function parameterization in weakly and moderately absorbing cases, while in strongly absorbing cases (mainly for β > 0.3) the accuracy deteriorates. Such strongly absorbing cases are, however, associated with small values of snow albedo and reflected radiances. The SSPs and the resulting snow albedo and reflected radiances for the OHC were compared with two previously used shape assumptions for snow grains, spheres and secondgeneration Koch fractals. The asymmetry parameter for the OHC is distinctly smaller than that for spheres but slightly higher than that for Koch fractals. Consistent with this, snow albedo for the OHC is generally substantially higher (slightly lower) than that for spheres (Koch fractals) for a given snow grain size r vp . Also for the distribution of reflected radiances, spheres differ more from the OHC than Koch fractals do.
The main limitation of the current work is that the SSP parameterization is based on a rather limited observational data set. The OHC was selected using scattering measurements at a single wavelength λ = 0.80 µm for only two cases with blowing snow. This raises several potential issues: -The choice of the OHC based on scattering measurements only implies that it most probably does not represent properly the actual distribution of snow grain shapes in blowing snow (or snow on ground). It also neglects the potential dependence of snow grain shapes on their size. Therefore, there is no guarantee that it represents the snow SSPs equally well at other wavelengths or for all snow grain sizes.
-Since absorption is very weak at λ = 0.80 µm, the observations do not constrain properly absorption by snow. Therefore, we cannot expect that our parameterization of β (Eq. 11) predicts precisely the actual values for snow. However, we do expect that it captures reasonably the systematic difference between non-spherical snow grains and spheres: in general β is larger for nonspherical particles.
-It is also possible that the snow grain shapes, and therefore the SSPs of snow on ground, might differ from those of blowing snow, and they might well vary from case to case, depending on how much metamorphosis the snow has experienced.
All these issues point to the need for validation of the derived parameterization against actual snow reflectance measurements in future work. In spite of the concerns mentioned above, it seems reasonable to assume that the OHC selected here provides a substantially better basis for representing the SSPs of snow than spheres do. Moreover, the parameterization equations provided in this paper are analytic and simple to use. A Fortran implementation of the snow SSP parameterizations is available at https://github.com/praisanen/snow_ssp.
To conclude, this paper describes a first-of-its-kind parameterization for representing the SSPs of snow in the solar spectral region. The parameterization is provided in hope that it will be useful, especially to those researchers that still use spherical particles for computing the radiative effects of snow. Nevertheless, it should definitely not be viewed as the "final solution" to the treatment of SSPs of snow. We hope that the present work will inspire the future development of snow SSP parameterizations based on more comprehensive data sets. Furthermore, at least in principle, it would be desirable to replace the current approach (where the shape distribution of snow grains is selected based on scattering measurements only) with an approach that more directly links the snow grain shapes to those actually observed. This would require, first, the parameterization of the size-shape distribution of snow grains based on observations and, second, the computation and parameterization of their SSPs. The main challenge in such an approach is the treatment of irregular grains, which are very common in snow.
www.the-cryosphere.net/9/1277/2015/ The Cryosphere, 9, 1277-1301, 2015  (Yang et al., 2013) β single-scattering co-albedo = 1 − single-scattering albedo δ Dirac's delta function θ zenith angle θ 0 zenith angle for incident radiation θ s scattering angle λ wavelength µ 0 cosine of zenith angle for incident radiation ν power-law index in the Legendre polynomial expansion of the correlation function of radius for Gaussian random spheres ξ non-dimensional absorption parameter (Eq. 8) σ relative SD of radius for Gaussian random spheres φ azimuth angle ω single-scattering albedo f truncated fraction of phase function in δ-M scaling (Wiscombe, 1977) g asymmetry parameter g 1 asymmetry parameter for the Henyey-Greenstein part in Eq. (21), defined by Eq. (24) g diff asymmetry parameter for diffraction (Eq. 19) g ray asymmetry parameter for the ray-tracing part (Eq. 23) m i imaginary part of refractive index m r real part of refractive index n max degree of truncation of the Legendre polynomial expansion of the correlation function of radius for Gaussian random spheres P projected area P 11 phase function P ref 11 reference phase function constructed from CLIMSLIP data (Eq. 2) P OHC 11 phase function for the optimized habit combination P HG Henyey-Greenstein phase function (Eqs. 18, 20) P diff parameterized phase function for diffraction (Eq. 17) P ray parameterized phase function for the ray tracing part (Eq. 21) P resid residual in the phase function parameterization (Eq. 25) P * resid residual in the phase function parameterization, truncated for δ-M scaling (Eqs. 26, 28) P n nth order Legendre polynomial Q ext extinction efficiency r vp volume-to-projected area equivalent radius t degree of distortion for Koch fractals V volume w 1 weight factor for the Henyey-Greenstein part in Eq. (21), defined by Eq. (22) w diff weight factor for the diffraction part in the parameterized phase function (Eqs. 14, 30), defined by Eq. (15) w ray weight factor for the ray tracing part in the parameterized phase function (Eqs. 14, 30), defined by Eq. (16) x size parameter x abs size parameter for absorption (Eq. 12) x p size parameter defined with respect to the projected area equivalent radius x vp size parameter defined with respect to the volume-to-projected area equivalent radius (Eq. 10)