Experimental and model-based investigation of the links between snow bidirectional reflectance and snow microstructure

Snow stands out from materials at the Earth’s surface owing to its unique optical properties. Snow optical properties are sensitive to the snow microstructure, triggering potent climate feedbacks. The impacts of snow microstructure on its optical properties such as reflectance are, to date, only partially understood. However, precise modelling of snow reflectance, particularly bidirectional reflectance, are required in many problems, e.g. to correctly process satellite data over snow-covered areas. This study presents a dataset that combines bidirectional reflectance measurements over 500–2500 nm and the X-ray tomography of the snow microstructure for three snow samples of two different morphological types. The dataset is used to evaluate the stereological approach from Malinka (2014) that relates snow optical properties to the chord length distribution in the snow microstructure. The mean chord length and specific surface area (SSA) retrieved with this approach from the albedo spectrum and those measured by the X-ray tomography are in excellent agreement. The analysis of the 3D images has shown that the random chords of the ice phase obey the gamma distribution with the shape parameter m taking the value approximately equal to or a little greater than 2. For weak and intermediate absorption (high and medium albedo), the simulated bidirectional reflectances reproduce the measured ones accurately but tend to slightly overestimate the anisotropy of the radiation. For such absorptions the use of the exponential law for the ice chord length distribution instead of the one measured with the X-ray tomography does not affect the simulated reflectance. In contrast, under high absorption (albedo of a few percent), snow microstructure and especially facet orientation at the surface play a significant role in the reflectance, particularly at oblique viewing and incidence.

1 Introduction

Snow optical properties are crucial to quantify the effect of snow cover on the Earth energy balance. They are also unique since snow is the most reflective material on the Earth surface. The subtle interplays between snow microstructure and snow optical properties are responsible for several climate feedbacks (e.g. Flanner et al., 2012). The dependencies of snow reflectance on snow microstructure have already been widely studied. It has been shown that the snow reflectance in the visible and near-infrared regions is primarily determined by the effective grain size that is defined as the ratio of the particle volume to its average projected area (Grenfell and Warren, 1999; Kokhanovsky and Zege, 2004) and is uniquely related to the specific surface area, hereafter SSA, defined as the ratio between the ice–air interface area and the mass of a snow sample (e.g. Flin et al., 2004; Domine et al., 2006; Gallet et al., 2009).

The effect of snow microstructure on the optical properties of snow is currently not fully understood. Up to now, many studies have focused on retrieving the single-scattering properties of individual ice crystals with “idealized” shapes (e.g. Xie et al., 2006; Picard et al., 2009; Liou et al., 2011; Räisänen et al., 2015; Dang et al., 2016) and on using these calculations to infer the effect of crystal shapes on snow optical properties. Several studies have already shown that the effect of shape is more pronounced on bidirectional reflectance (e.g. Dumont et al., 2010) and the vertical irradiance profile in the snowpack (e.g. Libois et al., 2013) than on hemispherical reflectance (albedo). Alternative approaches include running ray-tracing models directly on 3D images of the snow microstructure as done by Kaempfer et al. (2007) or at the intermediate level as done by Haussener et al. (2012) and Varsa et al. (2021). Similarly, Ishimoto et al. (2018) used X-ray tomography images of different snow types and a ray tracing model to compute the single-scattering properties of snow particles. They found that the modelled orientation-averaged scattering phase functions at two wavelengths (532 and 1242 nm) exhibit only a slight variation with the particle shapes.

Understanding and modelling the variations in snow directional reflectance with snow microstructure are essential to correctly interpret satellite data (Schaepman-Strub et al., 2006). Moreover, the sensitivity of snow directional reflectance to crystal shapes at least for high absorptive wavelengths (Xie et al., 2006; Dumont et al., 2010; Krol and Löwe, 2016a) makes snow directional reflectance a good candidate to provide an objective measurement of snow morphology. This characterization is often performed using the grain types defined in Fierz et al. (2009). Such an approach however depends on the observer and does not provide a variable that can be measured objectively and used directly in snowpack detailed modelling. Stanton et al. (2016) investigated the relationship between the bidirectional reflectance factor (BRF) and various snow crystal morphologies based on measurements in the 400–1300 nm range. They concluded that, as surface hoar grows, the snow surface becomes less and less Lambertian but that the geometry (illumination and viewing) at which the reflectance is maximum or minimum is difficult to predict. Horton and Jamieson (2017) used reflectance measurements and investigated the potentiality of normalized difference indices calculated from conical reflectance measurements at two wavelengths (860 and 1310 nm) to classify different crystal morphologies. They concluded that the bidirectional reflectance properties for different snow types must be investigated further.

The formalism developed by Kokhanovsky and Zege (2004) provides an analytical formulation linking the snow single-scattering properties and reflectance to the effective grain size and a shape parameter B, a.k.a. absorption enhancement parameter. This formulation has been used in the snow radiative scheme TARTES (Libois et al., 2013), enabling the retrieval of B values from concomitant measurements of irradiance profiles and reflectances on “homogeneous” snow layers. Though the obtained B values (from 1.4 to 1.8) significantly differ from the B value for spheres (1.25), no clear relationship has been established between grain types as defined in Fierz et al. (2009) and B (Libois et al., 2015). It must be underlined that this approach applies to low absorptions only, as B is the proportionality factor in the first (linear) term of the power expansion of the absorption coefficient of snow in the absorption coefficient of ice.

More recently, Malinka (2014) used the stochastic approach that considers a porous material as a random two-phase mixture and directly relates its optical properties to the chord length distribution (CLD) in the medium. This approach is not restricted to low absorptive wavelengths and directly relates the snow optical properties to the snow microstructure by means of the CLD; e.g. the shape parameter B for the random mixture arises in a natural way and equals n², where n is the refractive index of ice. For the random mixture with n=1.31, B=1.72, while measurements in real snow performed both in a laboratory and in situ gave values of B ranging from 1.4 to 1.8 (Libois et al., 2014). The approach has been successfully evaluated with respect to reflectance measurements over sea ice by Malinka et al. (2016).

In addition, Krol and Löwe (2016a) used X-ray tomography images to compare different metrics of snow microstructure. They experimentally demonstrated that the second moment of the CLD, μ₂, can be related to a curvature length and also theoretically to the absorption enhancement parameter B. They theoretically predicted, using the framework of Malinka (2014), that this microstructure metric strongly influences snow optical properties for high absorptive wavelengths. They also showed that the deviation of the CLD from the exponential law, which can be calculated using μ₂, varies with snow types. However, no measurement of snow optical properties was used in this study to evaluate the validity of their findings.

Figure 1 (a)–(c) Pictures of snow from S1–S3 taken with a microscope and (d) experimental set-up for S2 sampling. The picture shows the inner part of the temperature gradient box (Calonne et al., 2014a) and the metallic sample holder for the optical measurements.

To sum up, it has been shown that BRF, especially in high absorptive wavelengths, is more sensitive to snow morphology than snow albedo (bi-hemispherical reflectance). Yet no clear relationship has been either established theoretically or evaluated experimentally using optical measurements combined with an objective quantification of the snow microstructure. The objectives of the paper are thus to (i) describe one of the very few datasets that combined measurements of the bidirectional reflectance over the 500–2500 nm range and X-tomography characterization of the snow microstructure, (ii) evaluate the accuracy of the model of Malinka (2014) to simulate the snow BRF and its dependencies on snow microstructure parameters, and (iii) investigate the relationship between bidirectional reflectance and the snow microstructure beyond SSA using both simulation and measurement.

The first section provides a description of the snow samples, the measurement strategy, the optical model, the processing of the X-ray tomography images, and the optical data. The second section presents the results in terms of temporal variability of the snow microstructure, accuracy of the optical measurements, snow microstructure characterization, and model evaluation. Discussions and conclusions are detailed in the last section.

2 Materials and methods

2.1 Experimental set-up and sample description

The general idea of the experiment was to characterize both the snow microstructure and the BRF for the same macroscopic snow sample. The dataset consisted of three macroscopic snow samples: S1 analysed in March 2012 and S2 and S3 analysed in March 2013. S1 consisted of decomposing and fragmented particles/rounded grains (DF/RG) according to the classification of Fierz et al. (2009). S2 was composed of faceted crystals/depth hoar (FC/DH) obtained in a temperature gradient experiment. S3 was taken from the same temperature gradient experiment as S2 except that it was turned upside down so that the grain orientation was changed by 180^∘. Under temperature gradient, the facet formation is oriented toward the warmer side of the snow layer. For instance, when the temperature gradient is pointing downward, as is usually the case in nature, the facets tend to form on the downward surfaces while the upward surfaces stay more rounded (see e.g. Figs. 5 and 8 in Calonne et al., 2014a). As a consequence, by flipping S3, the faceted surfaces were oriented upward instead of downward in S2. This was done to investigate the effect of facet orientation on BRF. For each sample, the experiment involved the following steps (see Sect. 2.1.1–2.1.3 for more detail):

1. snow sample preparation,

2. snow microstructure characterization (manual measurements, casting for X-ray analysis) and preparation of a sample for BRF measurements,

3. BRF measurements,

4. snow microstructure characterization (manual measurements, casting for X-ray analysis).

Steps 2 and 4 were performed to characterize the microstructure just before and after the BRF measurements and to control the possible evolution of the microstructure during the BRF measurements.

2.1.1 Snow sample preparation

Figure 1 provides pictures of the snow from each sample. For S1, a 7 cm thick snow layer was collected on a 60 × 60 cm² styrodur plate after a snowfall close to the lab and stored for 3 weeks in isothermal conditions at −20 ^∘C. It then stayed 3 d at −10 ^∘C to reach the DF/RG state (Fig. 1a). The objective of this imposed isothermal metamorphism was to obtain a relatively recent snow sample, but with smooth and rounded shapes, that can be resolved at the pixel size we could access with the tomograph (between 6 and 12 µm). For S2 and S3 (Fig. 1b), fresh snow was collected in the field and was sieved into the temperature gradient box from Calonne et al. (2014a) (105 × 58 × 17.5 cm³). A vertical temperature gradient of ≈ 19.4 ^∘C m⁻¹ was applied inside the box with a mean temperature of −4 ^∘C. Such conditions produce simple structures of large and regular faceted crystals in a reasonable amount of time (Flin and Brzoska, 2008, and Calonne et al., 2014a). Before the measurements were made, S2 and S3 were stored under the described conditions for 16 and 18 d respectively.

2.1.2 Snow sampling and manual characterization

For the BRF measurements, the snow was sampled in a cylindrical sampler with no disturbance of the snow surface as in Dumont et al. (2010) (Fig. 1b). The diameter of the cylinder was 28 ± 1 cm. For S1 the snow thickness was 6.5 cm and for S2 and S3 it was 16.5 cm.

Before this sampling, the snow SSA and density were measured, the SSA with DUFISSS (DUal Frequency Integrating Sphere for Snow SSA measurement; Gallet et al., 2009) and ASSSAP (Alpine/Arctic Snow Specific Surface Area; Arnaud et al., 2011) and the density with manual weighing. The uncertainties on the SSA measured with DUFISSS and ASSSAP are in the range 10 %–15 % for SSA smaller than 60 m² kg⁻¹ (Gallet et al., 2009; Arnaud et al., 2011). For the density measured by manual weighing, the uncertainties are in the range of 1 % to 6 % (Proksch et al., 2016). Some small near-surface samples were taken for the X-ray analysis and casted using chloronaphthalene (chl) to block snow metamorphism (Flin et al., 2003; Calonne et al., 2014b). These measurements and sampling were reconducted after the optical measurements using the snow sampled in the big cylinder sampler.

2.1.3 X-ray tomography and BRF measurements

For all the samples, the X-ray tomography was performed at 7 µm resolution. Additionally, higher and lower resolutions (5 and 10 µm) were acquired for S1 and S2. The scanned samples were composed of three materials, namely ice, chloronaphthalene (chl), and some residual air bubbles due to incomplete impregnation of the sample (Flin et al., 2003; Hagenmuller et al., 2013). These three materials can be distinguished by their X-ray attenuation coefficient, i.e. their grayscale value I. Table 1 provides an overview of all the images taken for the three samples, and Fig. 2 shows a sub-sample of the 3D images obtained for each sample. The image name provides the sample name, the timing of the scan (B for “before the start of the BRF measurements”, A for “after the end of the BRF measurements”), the location of the sampling ((1) in the centre of the optical sample and (2) on the border of it), and the resolution in microns.

Table 1 Summary of all the X-ray tomography images acquired. “B” and “A” refer to “before” and “after” the optical measurement respectively.

Figure 2 Microstructure of the samples S1, S2, and S3 as revealed by X-ray tomography. These visualizations correspond to subsets from the 3D images S1_B_1_7m (a), S2_B_1_7m (b), and S3_B_1_7m (c). The scale bar is 1 mm.

The bidirectional reflectance was measured with a sensor field of view of 2.05^∘ using a set-up described in Dumont et al. (2010) and Brissaud et al. (2004) in a cold room at −10 ^∘C. The relative accuracy of the bidirectional reflectance distribution function (BRDF) measurements is estimated to be 1 % in Bonnefoy (2001). However, we do not believe that this accuracy is reached for high illumination angles (see Sect. 3.2).

Tables A1, A2, and A3 in Appendix A give an overview of the characteristics of the optical measurements for the three samples. The total duration of the optical measurements was 41 h for S1, 45 h for S2, and 94 h for S3.

2.2 X-ray tomography: image processing and analysis

2.2.1 Image processing

All grey level images were segmented using the following three-step semi-automatic method: (i) pre-processing of the image, including basic beam hardening and ring artefact corrections; (ii) detection of air bubbles and replacement of their levels by the mean grey value of 1-chloronaphthalene (see Flin et al., 2003, 2004, and “METHOD/Threshold-based segmentation/Air bubble detection” in Hagenmuller et al., 2013); and (iii) application of the energy-based binary segmentation method of Hagenmuller et al. (2013), with a resolution parameter r=1.2 voxel.

Once segmented, the obtained binary images can be described in terms of I(X), an indicator function of the ice phase such that

$$I\left(\mathit{X}\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if}\mathit{X}\text{lies in the ice phase}\\ \mathrm{0}& \text{if}\mathit{X}\text{lies in the air phase}\end{array},\right.$$

where $\mathit{X}=(x,y,z)$ is a position vector within the sample. Homogeneous cubic subsets of size n_max=700 voxels were then extracted from the original images for further analysis.

2.2.2 Chord length distribution

The chord length distribution (CLD), also called chord length probability density, is often used for the characterization of binary porous media. It is based on a microstructure description in terms of random chords, i.e. iso-phase line segments, whose lengths l are estimated by throwing virtual rays in random directions through the microstructure. The CLD p⁽ⁱ⁾(l) of phase i denotes the probability p⁽ⁱ⁾(l)dl of finding a random chord of length between l and l+dl in phase i (Torquato, 2002), thus giving us information on the thicknesses of the elements constituting the considered phase.

In the case of snow, which is known to be an anisotropic material with an orthotropic axis corresponding to the vertical (z) direction (see Calonne et al., 2011, 2012; Löwe et al., 2013; Calonne et al., 2014b, a; Wautier et al., 2015; Srivastava et al., 2016; Wautier et al., 2017), the CLD measured along a particular line depends on its direction. Assuming that the anisotropy is small, we consider the statistical characteristics of a sample in the three Cartesian directions. Namely, the chord lengths were obtained by scanning the segmented images with “rays” along the $\mathit{\{}\mathit{x},\mathit{y},\mathit{z}\mathit{\}}$ directions. The z axis is aligned with the vertical direction, while x, y are arbitrarily chosen in the plane perpendicular to the z axis. As the resolution of the X-ray images Δd is finite, the chord length l takes 700 discrete values from Δd to d_max=700Δd along every direction. The total number of chords of length l, $n_j^{\left(i\right)}\left(l\right)$, is then used to calculate the CLD of phase (i) in each X-ray image j:

$$\begin{array}{}\text{(1)}& p_j^{\left(i\right)}\left(l\right)={\displaystyle \frac{n_j^{\left(i\right)}\left(l\right)}{\mathrm{\Delta }{d}_j{\sum }_l{n}_j^{\left(i\right)}\left(l\right)}}.\end{array}$$

All chords that cross the file borders are inherently dismissed by this estimation method; hence, no hypothesis on the exact nature of the phase (air or ice) outside of the processed image is required. In what follows, the CLDs of the ice phase only are considered. The mean chord length a_j and the characteristic function (CF), L_j(α), of the ice phase in X-ray image j are obtained as

$$\begin{array}{}\text{(2)}& a_j=\langle l{\rangle }_j,\end{array}$$

$$\begin{array}{}\text{(3)}& L_j\left(\mathit{\alpha }\right)=\langle e^{-\mathit{\alpha }l}{\rangle }_j,\end{array}$$

where sign 〈〉_j denotes averaging in image j. The overall characteristics of the sample are the average over its X-ray images:

$$\begin{array}{}\text{(4)}& X={\displaystyle \frac{\mathrm{1}}{N}}\sum _j{X}_j,\end{array}$$

where X_j is any metric attributed to image j and N is a number of images of a sample (see Table 1 for the number of images per sample). The sample anisotropy is estimated as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{2}\langle l_z\rangle }{\langle l_x\rangle +\langle l_y\rangle }},\end{array}$$

where the subscript means the direction at which chords are taken.

2.2.3 Specific surface area (SSA)

The specific surface area was estimated using two different means: a stereological approach and a voxel projection method. Corresponding quantities calculated with these two methods are indicated with subscripts CLD and VP respectively.

The stereological approach is directly based on the computation of the mean chord length a of the ice phase obtained from the CLD analysis. Indeed, a is uniquely related to the SSA and to the ice volumetric mass ρ_ice by the following formula (see e.g. Torquato, 2002; Malinka, 2014):

$$\begin{array}{}\text{(6)}& {\text{SSA}}_{\text{CLD}}={\displaystyle \frac{\mathrm{4}}{a{\mathit{\rho }}_{\text{ice}}}}.\end{array}$$

The voxel projection estimation, denoted SSA_VP, corresponds to an approach developed by Flin et al. (2005, 2011). Based on an adaptive determination of the normal unit vector in each surface voxel, this method allows us to obtain the surface area of the whole object of interest. The SSA can then be obtained after dividing the computed interface area S_int by the associated volume. In the present implementation, a small improvement concerning the computation of the ice phase volume has been added. With the VP method being based on the area estimation of the interface located within a particular voxel (Flin et al., 2005), digitizing divides an image into voxels of two types: ice and air, with the interface voxels being systematically attributed to the ice phase. As a result, the total volume of ice voxels V_VP slightly overestimates the true volume of the ice phase V_ice by a quantity equal on average to a half volume of the interface voxels V_int. Therefore, SSA_VP was computed with the formula

$$\begin{array}{}\text{(7)}& {\text{SSA}}_{\text{VP}}={\displaystyle \frac{S_{\text{int}}}{{\mathit{\rho }}_{\text{ice}}\left(V_{\text{VP}}-\frac{V_{\text{int}}}{\mathrm{2}}\right)}},\end{array}$$

where S_int, V_VP, and V_int are computed from the X-ray images. This improved VP method has been validated against several data, including a series of calibrated spheres (see Appendix C).

2.3 Optical measurement analysis

As detailed in Dumont et al. (2010), to convert the spectrogonio-radiometer measurements into BRF values, we divide the measured radiance reflected from the snow surface by the radiance from a reference surface for which spectral albedo and BRF are known. For visible and near-IR wavelengths, the reference surface is a Spectralon^® panel. For wavelengths longer than 2440 nm, an infragold^® panel is used.

Let $R({\mathit{\theta }}_{\text{i}},{\mathit{\theta }}_{\text{v}},{\mathit{\varphi }}_{\text{i}},{\mathit{\varphi }}_{\text{v}},\mathit{\lambda })$ be the BRF of a sample under incident zenith angle θ_i, viewing zenith angle θ_v, incident azimuth ϕ_i, viewing azimuth ϕ_v, and wavelength λ. The reciprocity principle states that

$$\begin{array}{}\text{(8)}& R({\mathit{\theta }}_{\text{i}},{\mathit{\theta }}_{\text{v}},{\mathit{\varphi }}_{\text{v}}-{\mathit{\varphi }}_{\text{i}})=R({\mathit{\theta }}_{\text{v}},{\mathit{\theta }}_{\text{i}},{\mathit{\varphi }}_{\text{i}}-{\mathit{\varphi }}_{\text{v}}).\end{array}$$

In what follows, we use the anisotropy factor η calculated using the following equation to quantify the anisotropy of the reflectance of a snow sample:

$$\begin{array}{}\text{(9)}& \mathit{\eta }\left(\mathit{\lambda }\right)={\displaystyle \frac{R(\mathrm{30}{}^{\circ },\mathrm{70}{}^{\circ },\mathrm{180}{}^{\circ },\mathit{\lambda })-R(\mathrm{30}{}^{\circ },\mathrm{70}{}^{\circ },\mathrm{0}{}^{\circ },\mathit{\lambda })}{R(\mathrm{30}{}^{\circ },\mathrm{70}{}^{\circ },\mathrm{180}{}^{\circ },\mathit{\lambda })}}.\end{array}$$

Note that this parameter does not necessarily capture the position of the maximum and minimum reflectances, especially in the visible wavelengths as discussed in Stanton et al. (2016).

2.4 Optical modelling

The model of snow reflectance used in this investigation is described in detail in Malinka (2014) and Malinka et al. (2016). The concept notion used in this model to characterize the snow microstructure is the mean chord length a, which can be seen as the effective snow grain size (Eq. 6). The main quantity that characterizes the optical properties of a snow layer is its optical thickness, τ, which can be calculated as

$$\begin{array}{}\text{(10)}& \mathit{\tau }={\displaystyle \frac{\mathit{\beta }H}{a}},\end{array}$$

where β is the volume fraction of ice (1−β is the snow porosity) and H is the sample thickness.

Table 2 Sample microphysical properties calculated from the X-ray tomography images and retrieved from the spectral albedo. “B” and “A” refer to “before” and “after” the optical measurement respectively.

Further analysis showed that the albedo calculated with Eq. (10) was strongly overestimated in the green range of the spectrum. In particular, samples S2 and S3 with the geometrical thickness of 16.5 cm have the optical thickness τ calculated with Eq. (10) of greater than 200. According to Malinka et al. (2016) such an optical thickness produces in the green range an albedo of the order of $\frac{\mathit{\tau }}{\mathit{\tau }+\mathrm{4}}\approx \mathrm{0.98}$, while the measured quantities reliably show a value of about 0.9. This means that the snow samples contain some light-absorbing particles (Warren, 1982). These particles can be incorporated into the model. When their size is orders of magnitude smaller than that of snow grains (e.g. black carbon), scattering by these particles is negligible in comparison with scattering by snow grains. Thus, the effect of impurities can be modelled as an increase in the absorption coefficient of ice:

$$\begin{array}{}\text{(11)}& \mathit{\alpha }={\mathit{\alpha }}_{\text{ice}}+\mathit{\xi }{\mathit{\rho }}_{\text{ice}}{C}_{\text{m}},\end{array}$$

where α is the resulting effective absorption coefficient of ice, α_ice is the absorption coefficient of pure ice, ξ is the particle absorption cross section per its mass (mass absorption cross section), and C_m is the relative mass concentration of absorbing particles. The load and type of impurities were not measured in this experiment, so our choice of a pollutant was quite arbitrary. We assumed that the snow was polluted by black carbon with ξ=11.25 m² kg⁻¹ at 550 nm, the value recommended by Hadley and Kirchstetter (2012), and spectral dependence given by the inverse of wavelength λ⁻¹ (Bond and Bergstrom, 2006). The impurities internally mixed in the ice phase lead to an increase in the mass absorption efficiency (Flanner et al., 2012; He et al., 2018). The value of 11.25 m² kg⁻¹ is an intermediate value between fresh BC and internally mixed aged BC (Tuzet et al., 2019, 2020).

The model described in Malinka et al. (2016) uses the three parameters (the optical thickness τ, mean chord length a, and concentration of a contaminant with a predefined spectrum, C_m in our case) to characterize the snow reflectance in full. The more general model of light scattering in porous materials (Malinka, 2014) characterizes the medium microphysical properties via the chord length distribution (CLD). The model of snow in Malinka et al. (2016) assumes exponential CLDs. This model was shown to be reliable in terms of reproducing the measured albedo spectra with the simulated ones in the visible and near-IR ranges (from 350 to 1350 nm). However, in Malinka et al. (2016), there were no direct measurements of the snow microstructure characteristics to verify the model completely. In this work we verify the model by comparing the retrieved mean chord length with that measured directly by the X-ray tomography and investigate the bidirectional reflective properties in a rather wider spectral range (up to 3000 nm), addressing the effect of the µCT-measured CLD vs. the exponential one. To do so, we consider two configurations in the simulation.

• Assume that the CLD is exponential with the mean values calculated directly from the 3D images (Table 2, col. “Mean chord length, Aver.”). These simulations are labelled EXP hereinafter.

• Directly use the CLD calculated from the 3D images. The simulations are labelled µCT hereinafter.

For the simulations, we used the database of ice optical properties provided by Warren and Brandt (2008) except when another source is mentioned (Kou et al., 1993, and Grundy and Schmitt, 1998) in Sect. 3.4.

2.5 SSA retrieval from optical measurements

The optical model described in the previous section can be used to retrieve SSA from a measured albedo spectrum. First, the volume fraction of ice β is calculated for a sample:

$$\begin{array}{}\text{(12)}& \mathit{\beta }={\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\text{ice}}}},\end{array}$$

where ρ is the average sample density taken from the 3D X-ray images. Then the optical thickness τ is calculated by Eq. (10) with an arbitrary starting value of the mean chord length a. After that, the new values of a and C_m are found that provide the best fit, in the least-squares meaning, of the measured albedo spectrum and the simulated spectrum for a given τ calculated using the model from Malinka (2014). The new τ is then calculated by Eq. (10) with the new value of a, and the procedure is repeated until the resulting values no longer change in the first five digits. It usually takes two to three iterations. The SSA and a are related with Eq. (6).

3 Results

Here we first provide a comparison of the characterization of the snow microstructure obtained with the different methodologies (manual measurements, optical measurements, and X-ray tomography). In a second step, we evaluate the accuracy of the BRF measurements. The last two sections compare the BRF obtained from the model and from the measurements and evaluate the impact of the microstructure on the BRF.

3.1 Snow microstructure

3.1.1 SSA and density

The SSA was obtained by various methods using both optical measurements and X-ray tomography. The optical methods include DUFISSS, ASSSAP, and retrieval from the albedo spectrum (Sect. 2.3). The retrieved values are shown in Table 2 together with the values calculated from the 3D X-ray images. The measured and retrieved albedo spectra of samples S1 and S2 are shown in Fig. 3.

Figure 3 Measured and retrieved albedo spectra of samples S1 (a) and S2 (b). The albedo is retrieved following the method described in Sect. 2.5.

Figure 4 compares the SSA values obtained from the 3D images using either the chord length distribution calculation (SSA_CLD) or the voxel projection method (SSA_VP). It shows that the agreement between the two SSA values is good (R²=0.994, RMSE = 0.28 m² kg⁻¹). The lower SSA values are slightly overestimated via the CLD method, while the higher ones are slightly underestimated. This agreement is to be compared with Fig. 2 in Krol and Löwe (2016a).

Figure 4 SSA_CLD versus SSA_VP in m² kg⁻¹ computed for all 3D images described in Table 1. The triangle and circles markers indicate that the samples were taken before and after the optical scans respectively. Each marker corresponds to one X-ray image. Each macroscopic sample (S1–S3) is represented by a different colour. The linear fit is represented by the black dashed line.

Figure 5 sums up all the density and SSA measurements performed on all three samples.

Figure 5 Estimated density (first column) and SSA (second column) at various sampling depth in S1, S2, and S3. Values estimated from 3D images are represented by diamonds. Values estimated from measurements (manual weighting for density, DUFISSS, and ASSSAP for SSA) are represented by solid and dotted lines. The SSA retrieved from the spectral albedo is shown with vertical green lines. Values obtained before the optical scan are in blue and after in red. Note that the SSA obtained from 3D images represented in the figure is SSA_VP values. The snow surface is represented by the horizontal grey line. Numerical values of all measurements are reported in Tables D1 and D2.

All three samples exhibit some variability both horizontally and vertically. The average density estimated from the 3D images is slightly higher than the manually measured one for S1 and S3 and slightly lower for S2. The average SSA calculated from the X-ray tomography is systematically lower than that estimated from DUFISSS and ASSSAP measurements and is in perfect agreement with the SSA retrieved from the albedo spectra. These discrepancies might be inherent to the methodology but could also be linked to the variation in the properties inside the sample. Note that an X-ray tomography image has a very small size (of the order of several millimetres), the optical measurements with DUFISSS and ASSSAP use surfaces with the size of about 5 cm, and the spectral albedo characterizes a sample as a whole.

Figure 5 shows that S2 and S3 are denser and coarser – i.e. consist of larger grains – than S1 (decomposing and fragmented particles/rounded grains). The density after the BRF measurements is systematically higher than before, and the SSA is systematically smaller (except for one of the S3 measurements). This indicates that the snow microstructure as a whole has slightly evolved during the BRF measurements: snow has become denser and coarser (possibly because of sublimation in the cold room). This is to keep in mind when analysing the optical measurements. S1 exhibits larger changes in SSA than S2, and S2, in turn, demonstrates larger changes than S3. This is compatible with the fact that snow with lower density and higher SSA is subject to more rapid evolution (e.g. Flin et al., 2004; Carmagnola et al., 2014; Schleef et al., 2014). Also, the retrieval has shown that the snow taken in the mountains (S2 and S3) was much more pure: the soot concentration in it was 5 times less than that in the snow taken in the urban area (S1) (Table 2).

Figure 6 (a) CLDs in the three samples (circles and dotted curves) and their approximations with the gamma distribution with the parameter m (solid lines). (b) Same as the left panel but using the log scale. All the images from Table 1 have been used.

3.1.2 CLD analysis

Figure 6 compares the CLD of the three samples estimated from all the images. S1’s CLD is narrower than S2's and S3s' CLD; however, all of them exhibit similar features. They have almost exponential tales for large l and approach zero at l=0. The latter fact means the ice–air interface does not have edges at the image resolution for the three snow samples, in contrast, e.g. to a set of random polyhedra, which has a Markov property, i.e. the exponential CLD. Distributions of this type can be represented by the gamma distribution:

$$\begin{array}{}\text{(13)}& p\left(l\right)={\displaystyle \frac{l^{m-\mathrm{1}}}{\mathrm{\Gamma }\left(m\right)}}{\left({\displaystyle \frac{m}{a}}\right)}^m\exp \left(-{\displaystyle \frac{m}{a}}l\right),\end{array}$$

where a is an average, Γ(m) is the gamma function, and m is a shape parameter. The exponential distribution is a particular case of the gamma distribution with m=1. Thus, the shape parameter m indicates the deviation of a distribution of this type from the exponential one. Fitted gamma distributions are also indicated in Fig. 6 (see below for the calculation of the value of m).

The CLD determines optical properties of a mixture not explicitly but via its Laplace transform (Malinka, 2014), which by definition is the characteristic function of the distribution, L(α) (see Eq. 3).

$$\begin{array}{}\text{(14)}& L\left(\mathit{\alpha }\right)=\langle e^{-\mathit{\alpha }l}\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-\mathit{\alpha }l}p\left(l\right)dl\end{array}$$

If the argument α is equal to the substance absorption coefficient, then the characteristic function L(α) describes the process of photon absorption while travelling along random chords within an absorbing material, in this case ice (see Malinka, 2014). The Laplace transform of the gamma distribution (Eq. 13) is

$$\begin{array}{}\text{(15)}& L\left(\mathit{\alpha }\right)={\left(\mathrm{1}+{\displaystyle \frac{\mathit{\alpha }a}{m}}\right)}^{-m}.\end{array}$$

Figure 7 µCT-measured characteristic functions for S1 and S3 and their approximations with the Laplace transform of the gamma distribution with different shape parameters m. All the images have been used for S1 and S3. S2 is not shown because of the overlap with S3.

Figure 7 shows the characteristic functions calculated directly from the 3D images of S1 and S3 using Eq. 3 and their approximations with the Laplace transforms (Eq. 15) of the gamma distribution with m calculated by least-squares fitting. The Laplace transform of the exponential distribution (m=1) is shown for comparison. The characteristic function of S2 is not shown because of the overlap with that of S3.

The characteristic functions calculated directly from the X-ray images at large values of argument αa could be distorted by the image discretization. Their values are reliable if the error of the argument αa due to the discretization is small. If the image resolution is Δd, then the condition Δd≪1 should hold true. If the resolution is 7 µm and the mean chord is about 200 µm, we have

$$\begin{array}{}\text{(16)}& \mathit{\alpha }\mathrm{\Delta }d=\mathit{\alpha }a{\displaystyle \frac{\mathrm{\Delta }d}{a}}\approx {\displaystyle \frac{\mathit{\alpha }a}{\mathrm{30}}},\end{array}$$

i.e. αΔd<0.1 for αa<3, so the µCT-measured characteristic functions presented in Fig. 7 are not affected by the discretization error.

The characteristic function in Eq. (3) has the power series expansion

$$\begin{array}{}\text{(17)}& \begin{array}{rl}L\left(\mathit{\alpha }\right)& =\langle \mathrm{1}-\mathit{\alpha }l+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{\alpha }l{)}^{\mathrm{2}}+\mathcal{O}({\mathit{\alpha }}^{\mathrm{3}}{l}^{\mathrm{3}})\rangle \\ & =\mathrm{1}-\mathit{\alpha }a+{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{2}}}+\mathcal{O}\left({\mathit{\alpha }}^{\mathrm{3}}{a}^{\mathrm{3}}\right)\end{array},\end{array}$$

where μ₂ is the second moment of the CLD.

The exponential CLD has the Laplace transform from Eq. (15) with m=1, hence the expansion

$$\begin{array}{}\text{(18)}& L\left(\mathit{\alpha }\right)=\mathrm{1}-\mathit{\alpha }a+{\mathit{\alpha }}^{\mathrm{2}}{a}^{\mathrm{2}}+\mathcal{O}\left({\mathit{\alpha }}^{\mathrm{3}}\right).\end{array}$$

Comparing the second-order term in Eq. (18) with that in the general expansion in Eq. (17), we find that the deviation of the CLD from the experimental law, in general, can be characterized by the value $\frac{{\mathit{\mu }}_{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}}$ (Krol and Löwe, 2016a) or its deviation from unity, δ:

$$\begin{array}{}\text{(19)}& \mathit{\delta }=\mathrm{1}-{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}}}.\end{array}$$

For the gamma distribution this value is expressed via the shape parameter as follows:

$$\begin{array}{}\text{(20)}& \mathit{\delta }={\displaystyle \frac{m-\mathrm{1}}{\mathrm{2}m}}.\end{array}$$

The deviations from the exponential law δ, calculated from m values with Eq. (20) and directly from the CLD with Eq. (19), are shown in Table 3. The values obtained in this experiment are slightly higher but consistent with the results of the analysis of Krol and Löwe (2016a).

Table 3 Values of m and δ averaged across all the images.

The deviation of the CLD from the exponential law matters, i.e. affects the optical properties, if the difference between Eqs. (17) and (18) is not negligible, i.e. ${\mathit{\alpha }}^{\mathrm{2}}{a}^{\mathrm{2}}-\frac{{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{2}}}{\mathrm{2}}={\mathit{\alpha }}^{\mathrm{2}}{a}^{\mathrm{2}}\mathit{\delta }\sim \mathrm{1}$.

Thus, the absorption at which the deviation matters can be estimated as

$$\begin{array}{}\text{(21)}& \mathit{\alpha }\sim {\displaystyle \frac{\mathrm{1}}{a\sqrt{\mathit{\delta }}}}.\end{array}$$

Provided a∼200 µm and δ∼0.25, we get α∼10⁴ m⁻¹. So, we can expect that the true shape of the CLD will impact the optical properties only at very high absorptions. The absorption coefficient of ice in the range 500–2500 nm approaches such values at 1500, 2000, and 2500 nm (see below). At these absorptions, for which αa∼2 , we can expect that the light penetration depth will not exceed the mean chord length a, treated as the effective snow grain size, so all the incident light will be absorbed within the skin layer – a layer at the surface with a thickness of “one grain”. The reflectance will be completely determined with Fresnel reflection by crystal facets and refraction by fine grains in this skin layer. Under these conditions, the fine grains, and therefore the CLD behaviour at small lengths, will play a more important role. The facet orientation at the surface can differ from that in deeper layers as well.

Figure 8 Measured reflectances of the three samples (each panel) obtained for optically equivalent geometries.

Figure 9 Measured (circles) and modelled (curves) BRF for the three samples (S1: a, S2: b, S3: c) for ${\mathit{\theta }}_{\text{i}}=\mathrm{30}{}^{\circ }$, ${\mathit{\theta }}_{\text{v}}=\mathrm{60}{}^{\circ }$, and ${\mathit{\varphi }}_{\text{v}}=\mathrm{180}{}^{\circ }$. The solid and dashed curves correspond to the simulated reflectance using the snow microstructure measured before and after the optical measurements respectively. All the simulations were performed with the exponential CLD. The differences between the measured and simulated reflectances are reported on the right y axis with crosses and dots. The 3D images at 7 µm resolution, i.e. three images per sample, are used for the simulations.

3.2 Optical measurement accuracy: reciprocity principle

Figure 8 compares the reflectances obtained for optically equivalent geometries according to Eq. (8). For S1, Eq. (8) is really well verified. This is not the case for S2 and S3, especially in the visible range. The SSA temporal evolution during the BRF measurements (see Fig. 5 in Sect. 3.1.1) cannot explain these discrepancies, since it would result in a decrease in reflectance in the measurements taken later, so that the red curves would be lower than the blue ones, but the contrary is observed. A possible explanation lies in the fact that the illumination pattern is larger in the sample holder for large incidence angles (e.g. 60^∘). Since the snow does not perfectly fill the sample holder, direct reflection of light coming in on the side of the sample holder to snow cannot be avoided, acting as an additional source of light and contributing to the reflected signal. This hypothesis is reinforced when comparing other equivalent geometries (θ_i=0^∘, θ_v=30^∘) showing almost perfect agreement between the two reflectances. This indicates that reflectances measured for θ_i=60^∘ should be analysed accounting for the fact that they might be overestimated because of re-illumination, especially for low absorptive wavelengths.

This shows that the relative accuracy of the measurements estimated at 1 % in Brissaud et al. (2004) is not reached for high illumination angles. In the following, the interpretation of the BRDF measurements is mostly restricted to wavelengths larger than 800 nm and – except in Sect. 3.4 – to an illumination angle of 30^∘ to minimize this effect. For wavelengths larger than 800 nm, the effect of the limited thickness of the sample and of the presence of light-absorbing impurities (Warren, 1982) is also limited.

3.3 Model vs. measurements

3.3.1 Spectral variations

Figure 9 compares the measured and simulated spectral reflectances for all three samples. The results show that the simulated reflectances generally agree well (absolute difference less than 0.02) except for several wavelength ranges that exhibit systematic bias for all three samples. There is a little BRF overestimation in the visible and near-IR ranges up to 1400 nm. Nevertheless, the albedo does not demonstrate any overestimation in this range of wavelengths. The BRF is overestimated mainly at oblique incidence and viewing angles and only in the visible range (see also Figs. 13, B3, and B4). We attribute the overestimation either to shortcomings of the measurements due to geometry (see above) or to the drawbacks of the model in the description of the angular dependence of the bidirectional reflectance, which may not take into account several factors, such as dense packing or surface roughness. There are some peaks of the difference – 1000, 1200, and 1450 nm – where the reflectance, following ice absorption, has steep changes. However, the most pronounced differences are in the ranges 1400–1900 and 2150–2300 nm. These differences are especially pronounced for S1. We relate these differences to the difficulties of the measurements of ice absorption in these ranges with high accuracy.

To examine this effect, we simulate the snow reflectance in the most pronounced difference range 1400–1900 nm using the values of the complex refractive index of ice taken from different databases. Figure 10 compares the measured and the simulated reflectances in this range for S1 and S2, using several databases of the complex refractive index of ice: Warren (1984), Kou et al. (1993), Grundy and Schmitt (1998), and Warren and Brandt (2008). For S1 both the incidence and observation angles are 60^∘. For S2 the observation angle is 60^∘, but the incident angle is taken as 30^∘, because of the problems discussed in Sect. 3.2. As seen from Fig. 9, the data from Grundy and Schmitt (1998) provide the simulated values closest to the measured ones and will be used in further simulations for this spectral region. However, as far as we understand, the additional careful measurements of the complex refractive index of ice in the IR range are still needed.

Figure 10 Measured (circles) and simulated (EXP, lines) reflectances in the principal plane for S1 (a) and S2 (b) in the 1400–1900 nm range. The incident zenith angles are 60^∘ for S1 and 30^∘ for S2. The viewing angle is 60^∘ zenith and 180^∘ azimuth. The BRF has been simulated using different ice refractive index databases: Warren and Brandt (2008) (light blue), Warren (1984) (dark blue), Grundy and Schmitt (1998) (red), and Kou et al. (1993) (orange).

Figure 11 Measured (circles) and simulated (EXP) (curves) BRF of the three samples at vertical incidence at 800 nm.

Figure 12 Measured (circles) and simulated (EXP, curves) BRF for S1 (a) and S2 (b) at vertical incidence at 1300 nm.

Figure 13 BRF measured (circles) and simulated with the exponential (curves) and µCT-retrieved (dashes) CLDs for S1 at 700, 1030, 1300, and 1500 nm. The incident zenith angle is 30^∘. The different colours correspond to the azimuth. The µCT-retrieved CLDs are averaged across all of the S1 X-ray images.

3.3.2 Angular distribution

Figures 11 and 12 compare the simulated and measured BRF at vertical incidence for 800 and 1300 nm. The simulated and measured data agree well, except for S1 at the observation directions close to vertical. Figure 13 shows the simulated and measured BRF at 30^∘ incidence for a wide range of the observation angles. Except at 1500 nm, the simulated BRF generally agrees with measured values. However, for high viewing zenith angle, typically higher than 45^∘, the simulated anisotropy is stronger than the measured one. These conclusions hold for S2 and S3 (see Figs. B3 and B4 in the appendices). At 1500 nm, the simulated BRF is lower than the measured values, as already reported in Fig. 9. The discrepancies between the simulated and measured BRF are generally greater than the spread in the simulation due to the variability of the snow microstructure (different images).

Figures 14, B1, and B2 (see Appendices) compare the measured and simulated angular distribution of reflectance for three wavelengths (1320, 800, and 1500 nm), for an incident zenith angle of 30^∘ and for the two configurations of simulations (EXP and µCT).

Figure 14 is obtained for 1320 m (wavelength with medium absorption). It shows that (i) the angular pattern and reflectance magnitude are reasonably reproduced by the model and (ii) there are no detectable differences between the EXP and the µCT simulations. It also shows that the higher reflectance values are obtained for S1 and the lower for S3, which is consistent with the SSA of the samples. The third line of Fig. 14 demonstrates that for the three samples, the simulated reflectance is overestimated in the forward scattering direction and underestimated in the backward direction, so at 1320 nm, the model seems to overestimate the anisotropy of the reflected radiation. The conclusions obtained at 1320 nm also stand at 800 nm (low absorption, Fig. B1). The possible causes of this general overestimation of the anisotropy of the angular pattern are discussed in Sect. 4. However at 1500 nm (high absorption, Fig. B2), the reflectance is underestimated almost twice for all angles and all samples in the simulations. The reflectances for this wavelength are very low and can also be affected by experimental uncertainties. Stronger differences in reflectances at 1500 nm are noticeable between the samples and also between the EXP and µCT simulations.

Figure 14 Measured (first line) and simulated (second line) reflectances for the three samples (three columns) at 1320 nm. The incident zenith angle is 30^∘. On the polar plot, the polar angle represents the azimuth, and the radius is proportional to the viewing angles. The dotted circles represent viewing angles of 30, 60, and 70^∘. The two lower lines represent the difference between the simulated and measured reflectance for the exponential and µCT-retrieved CLDs. The images used for the simulations are S1_E_1_7m, S2_B_1_7m, and S3_B_1_7m. .

Figure 15 Measured (crosses) and simulated (dotted curves for µCT and solid curves for EXP) values of parameter η (Eq. 9) obtained for the three samples. All the images were used for the simulations.

Figure 15 compares the anisotropy of the reflectance measured and simulated quantified by the anisotropy parameter η (Eq. 9). It shows that η is higher in the simulations than in the measurements, except for highly absorptive wavelengths (around 1500 and 2000 nm), where the differences between the samples are higher in the measurements than in the model. Also, the EXP and µCT simulations only differ for these highly absorptive wavelengths. This general overestimation of the anisotropy in the simulations is also visible in Figs. 9, 13, 14, B1, B2, B3, and B4.

Figure 16 Measured (markers) and simulated (lines) reflectance in the principal plane (forward direction) for an incident zenith angle of 60^∘ at 1280, 1420, 1500, and 2000 nm. Plain lines correspond to EXP and dotted lines to µCT simulations. All images have been used as inputs for the simulations. The red and orange lines correspond to simulation obtained for S1 with refractive index from Grundy and Schmitt (1998) and Kou et al. (1993). For all the other simulations we used refractive index from Warren and Brandt (2008) .

3.4 Results at 60^∘ illumination angle

Figure 16 compares the measured and simulated reflectance in the principal plane (forward direction) for an incident angle of 60^∘ at several wavelengths. At high viewing angles, for all wavelengths and despite the lower SSA, S3 exhibits higher reflectance than S2 and S1. The higher the absorption, the more pronounced this effect, so it might be related to the facet orientation at the surface. At 1280 nm, the simulations and the measurements agree well, with a slight overestimation of the reflectance at high viewing angles. For the other wavelengths, the model generally underestimates reflectance. The use of alternative ice refractive index values (Grundy and Schmitt, 1998) improves the simulations at 1420 and 1500 nm but is not sufficient to reconcile the measurements and the simulations for high viewing angles at 1500 and 2000 nm (see Sect. 4 for additional discussion). The difference between the µCT and EXP simulations is only noticeable for high viewing angles and high absorption. Under these conditions, accounting for the CLD shape retrieved from µCT leads to a decrease in the calculated reflectance in comparison with the simulations with exponential CLD.

4 Discussions and conclusions

This study presents a dataset that combines the snow bidirectional reflectance over the 500–2500 nm range with different illumination geometries and 3D images of snow using X-ray tomography, which allows analysis of the snow microstructure, e.g. its SSA and density. The SSA is calculated using two different methods: the voxel projection method (VP) and the method based on the ice chord length distribution (CLD). The comparison between the two SSA values is in excellent agreement (R²=0.994). The SSA and mean chord values computed from 3D images are generally lower than those obtained via DUFISSS and ASSSAP optical methods and in excellent agreement with those retrieved via spectral albedo. The comparison between density values obtained via the 3D images and via manual measurements exhibits slightly higher density for 3D images, which might relate to the heterogeneity of the samples and to the segmentation process applied to the images (see Fig. 12 in Hagenmuller et al., 2013). The analysis of the 3D images has shown that the CLDs approach zero at l=0 both in the DF/RG snow and in the snow after long evolution in temperature gradient conditions. This means that the air–ice interface is smooth and has no sharp edges even though faceted crystals are present, at least for the images considered in this study. A possible reason for this is the continued curvature-driven metamorphism of snow, which already begins during a snowfall (Flin et al., 2003; Brzoska et al., 2008; Krol and Löwe, 2016b).

The analysis of the characteristic functions of random chords in the snow phase calculated directly from the 3D images has shown that the random chords obey the gamma distribution with the shape parameter m, equal to 2.40, 1.95, and 1.94 for S1, S2, and S3 respectively. The deviation of the distribution from the exponential one is the lowest for the more faceted crystals (sample 3) and the highest for DF/RG (sample 1). Its values, however, largely vary for the same macroscopic sample probably due to both the temporal evolution and the snow heterogeneity. The deviation values are slightly higher but in the same range as obtained by Krol and Löwe (2016a).

The comparison between the simulated and measured reflectance under a specific geometry (${\mathit{\theta }}_{\text{i}}=\mathrm{30}{}^{\circ }$ or 60^∘, ${\mathit{\theta }}_{\text{v}}=\mathrm{60}{}^{\circ }$ and ${\mathit{\varphi }}_{\text{v}}=\mathrm{180}{}^{\circ }$) shown in Fig. 9 demonstrated that the simulated values generally agree well (absolute difference lower than 0.03) over the whole spectrum 500–2500 nm and that for this geometry the impact of the CLD is small. Systematic differences are found in several wavelength ranges. Such discrepancies have already been reported in several studies (e.g. Carmagnola et al., 2013) and attributed to uncertainties in the values of the ice complex refractive index from Warren and Brandt (2008). The use of alternative databases, especially the one from Grundy and Schmitt (1998) noticeably improves the agreement in the range 1400–1500 nm as shown by Fig. 10. A more extensive comparison of the simulated and measured reflectance angular distribution at 30^∘ illumination angle shows that the BRF predicted by the model is in good agreement with the measured one. Spectral changes in the angular distribution are well reproduced. The impact of the CLD shape on the simulation is only detectable for absorptive wavelengths (1500 nm) and high viewing zenith angles. The anisotropy of the reflectance, quantified by the relative difference between the reflectances measured at a 70^∘ zenith angle in the forward and backward directions in the principal plane, is systematically overestimated in the simulations. This effect was reported by many authors and commonly attributed to the surface roughness (Warren et al., 1998; Painter and Dozier, 2004; Hudson et al., 2006; Jin et al., 2008; Carlsen et al., 2020), which was not accounted for in our simulations and generally leads to less anisotropic angular patterns. However, the opposite effect is observed in the absorption bands around 1500 and 2000 nm, where the larger discrepancies are found between the samples. For these wavelengths the model underestimates the intensity of forward scattering in the principal plane. Despite the lower SSA, sample 3 exhibits the higher reflectance at 70 and 75^∘ zenith viewing angles. This might be related to the preferred upward orientation of the facets at the surface of this sample.

To sum up, the results exhibit two different trends for small and medium/long optical paths in the snow.

• For long/medium optical paths, the model predicts no impact of the CLD shape on the reflectance and is in good agreement with the measurement. The anisotropy is slightly overestimated in the simulations. A possible explanation can be related to the surface roughness, which is not accounted for in the simulations and generally leads to less anisotropic patterns (Hudson et al., 2006). Differences between the reflectance of the three samples in the IR are well reproduced using the exponential CLD.

• The picture is different for the small optical paths, which refer to the high absorptive wavelengths and to oblique illumination and/or viewing. The model predicts an impact of the CLD. The use of the µCT-retrieved CLD instead of the exponential distribution leads to a decrease in the intensity of the forward scattering peak in the principal plane. The higher forward scattering intensities are observed for the more faceted crystals (S3) with facets upward. S3 also depicts the CLD closest to the exponential one. So, this behaviour is consistent with the model prediction. Several hypotheses can be drawn to explain the model–measurements discrepancies. The first hypothesis is related to the imaginary part of the ice refractive index, to which the reflectance, including albedo, is very sensitive, but it is probably not sufficient to explain the different behaviours of the bidirectional reflectance distribution shown in Figs. 15 and 16. The second possible reason is the anisotropy of snow, which was noted by several authors (Calonne et al., 2012; Löwe et al., 2013) and for our samples ranged from 0.9 to 1.07, being calculated with Eq. (5). The last but most probable reason is the “skin effect”. For small optical paths, the light penetration depth is really small and probably less than a millimetre (e.g. Mary et al., 2013), so the role of the skin layer, where the grains could be finer and have a preferable facet orientation, is crucial. This is also corroborated by the fact that S2 and S3 have only a small difference in SSA, measured in deeper layers, but a strongly different behaviour in reflectance.

The effect of crystal shape and SSA on snow BRF is intricate (e.g. Dumont et al., 2010; Stanton et al., 2016). The dataset presented in this study is limited to only three snow samples, so no statistically robust conclusions on the effect of crystal shape on reflectance can be drawn. The measurements, however, show that the faceted crystals exhibit more anisotropic reflectance than fragmented particles/rounded grains as observed in Stanton et al. (2016) and that the anisotropic behaviour is further reinforced if facets are orientated upward. More quantitative relationships between crystal shape and the BRF would require a much larger number of snow samples analysed. One way of addressing this problem is to take advantage of ray tracing models such as those developed by Kaempfer et al. (2007), Picard et al. (2009), and Petrasch et al. (2007), which can be run directly on 3D images of the microstructure.

To conclude, this unique dataset combining X-ray tomography imaging of snow microstructure and high-accuracy measurements of snow BRF was used to demonstrate the following.

• Faceted crystals exhibit a more anisotropic reflectance than fragmented particles.

• The Malinka et al. (2016) model can be used to accurately simulate the snow BRF using SSA as input. The model, however, slightly overestimates the snow reflectance anisotropy. Even so, the mean chord length and SSA retrieved from the albedo spectrum and those measured by the X-ray tomography are in excellent agreement. As far as we know, such a successful comparison of the mean chord and SSA of snow retrieved from optical and µCT measurements has been obtained for the first time.

• Other characteristics of snow microstructure besides SSA, e.g. the CLD shape, impact the angular reflectance of snow for high ice absorption and oblique viewing and illumination.

Appendix A: Overview of the optical measurements

Table A1 Geometry and spectral range of S1 optical measurements.

Table A2 Geometry and spectral range of S2 optical measurements.

Table A3 Geometry and spectral range of S3 optical measurements.

Appendix B: Model vs. measurements

Figure B1 Measured (first line) and simulated (second line) reflectances for the three samples (three columns) at 800 nm. The incident zenith angle is 30^∘. On the polar plot, the polar angle represents the azimuth, and the radius is proportional to the viewing angles. The dotted circles represent viewing angles of 30, 60, and 70^∘. The two lower lines represent the difference between the simulated and measured reflectance for the exponential and μCT-retrieved CLDs. The images used for the simulations are S1_E_1_7m, S2_B_1_7m, and S3_B_1_7m.

Figure B2 Measured (first line) and simulated (second line) reflectances for the three samples (three columns) at 1500 nm. The incident zenith angle is 30^∘. On the polar plot, the polar angle represents the azimuth, and the radius is proportional to the viewing angles. The dotted circles represent viewing angles of 30, 60, and 70^∘. The two lower lines represent the difference between the simulated and measured reflectance for the exponential and μCT-retrieved CLDs. The images used for the simulations are S1_E_1_7m, S2_B_1_7m, and S3_B_1_7m.

Figure B3 BRF measured (circles) and simulated with the exponential (curves) and μCT-retrieved (dashes) CLDs for S1 at 700, 1030, 1300, and 1500 nm. The incident zenith angle is 30^∘. The different colours correspond to the azimuth. The μCT-retrieved CLDs are averaged across all of the S2 X-ray images.

Figure B4 BRF measured (circles) and simulated with the exponential (curves) and μCT-retrieved (dashes) CLDs for S1 at 700, 1030, 1300, and 1500 nm. The incident zenith angle is 30^∘. The different colours correspond to the azimuth. The μCT-retrieved CLDs are averaged across all of the S3 X-ray images.

Appendix C: Evaluation of SSA on calibrated spheres using different methods

Figure C1 SSA estimated from different methods for spheres and their relative error compared to the theoretical value SSA_th. SSA_CLD and SSA_VP are defined in Eqs. (6) and (7) respectively. SSA_{VP_2011} is described in Flin et al. (2011). Sphere radius ranges from R=1 to R=140 voxel (1 voxel was arbitrarily chosen to correspond to 10 µm) – overview in log scale. Note that here, only individual spheres (perfectly isotropic objects) are considered, which favours the SSA_CLD estimation.

Appendix D: Density and SSA manual measurements

Table D1 Manual measurements of density. The heights of the upper and lower boundaries of the sample used for density measurements are indicated in the table.

Table D2 Manual measurements of SSA. The heights of the upper and lower boundaries of the sample used for SSA measurements are indicated in the table.