# Unraveling the optical shape of snow

Alvaro Robledano, Ghislain Picard, Marie Dumont, Frédéric Flin, Laurent Arnaud, Quentin Libois

<figure><img src="image_1.png" /></figure>

To cite this version:

Alvaro Robledano, Ghislain Picard, Marie Dumont, Frédéric Flin, Laurent Arnaud, et al. Unraveling the optical shape of snow. Nature Communications, 2023, 14 (1), pp.3955. 10.1038/s41467-023-39671-3. hal-04245588

HAL Id: hal-04245588

https://hal.science/hal-04245588v1

Submitted on 20 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Unraveling the optical shape of snow

Received: 28 February 2023

Accepted: 21 June 2023

Published online: 07 July 2023

<figure><img src="image_4.png" /></figure>

Check for updates

Alvaro Robledano Ⓔ<fnref n="1" />,²✉, Ghislain Picard<fnref n="1" />, Marie Dumont Ⓔ², Frédéric Flin Ⓔ², Laurent Arnaud<fnref n="1" /> & Quentin Libois Ⓔ³

The reflection of sunlight off the snow is a major driver of the Earth's climate. This reflection is governed by the shape and arrangement of ice crystals at the micrometer scale, called snow microstructure. However, snow optical models overlook the complexity of this microstructure by using simple shapes, and mainly spheres. The use of these various shapes leads to large uncertainties in climate modeling, which could reach 1.2 K in global air temperature. Here, we accurately simulate light propagation in three-dimensional images of natural snow at the micrometer scale, revealing the optical shape of snow. This optical shape is neither spherical nor close to the other idealized shapes commonly used in models. Instead, it more closely approximates a collection of convex particles without symmetry. Besides providing a more realistic representation of snow in the visible and near-infrared spectral region (400 to 1400 nm), this breakthrough can be directly used in climate models, reducing by 3 the uncertainties in global air temperature related to the optical shape of snow.

Ice crystals formed in the atmosphere show a large variety of sophisticated and, often, near-perfect geometric shapes<fnref n="1" /> (Fig. 1a). The interaction of sunlight with such crystals sometimes results in well-known optical phenomena, called halos², whose nature is directly related to the shape of the crystals³.

Once the ice crystals are deposited on the ground and form the snow cover, they establish bonds with each other and their shape drastically changes. As a consequence, the original atmospheric crystals are rarely distinguishable after a few hours or days⁴,⁵, making snow cover a two-phase porous material made of ice and air, rather than a collection of individual particles. Snow is continuously evolving due to the thermodynamical non-equilibrium between ice and interstitial air that leads to recrystallisation of the water molecules, a process called metamorphism. This process leads to a variety of snow morphologies at the micrometer scale, called microstructure hereinafter, which are less regular, less symmetrical and more diverse than the original crystals (Fig. 1b–d). Hence the interaction of sunlight with snow on the ground is more complicated to model than with ice clouds in the atmosphere, and is not yet fully understood. The resulting uncertainties significantly alter the estimation of the solar radiation reflected by snow and in turn the Earth's radiative budget⁶,⁷. Recent climate simulations thus show that simply changing the snow microstructure can modify the global annual-mean 2 m air temperature by nearly 1.2 K⁸.

Early studies have tried to relate the snow optical properties to the snow microstructure, relying on strong approximations such as considering snow as a collection of disconnected and independent ice spheres⁹. Under this hypothesis, the Mie theory provides the optical properties for a given sphere size. Later on, the relative success of the equivalent-sphere concept, i.e., the representation of snow as a collection of spheres with the same volume-to-area ratio, to simulate snow albedo, led to the spread of this concept in most snow optical models¹⁰⁻¹². This representation provided a clear relationship between the snow albedo and the sphere size¹⁰,¹³, called the “optical diameter”. However, posterior studies highlighted several caveats of the spherical assumption¹⁴⁻¹⁶, particularly in determining how the snow reflects light in different directions and how deep the light penetrates¹⁷⁻²⁰. This uncertain representation has implications for the interpretation of satellite data²¹,²², for snow photochemistry²³⁻²⁵, for light transmission through snow over sea ice²⁶, and more importantly for the surface energy budget of snow-covered regions²⁶⁻²⁸. This raises the need to go beyond the unrealistic but still widely used spherical representation of snow.

Many alternative shapes have been used to describe the snow microstructure, such as fractals²⁷, Voronoi tessellations²⁸, cubes, hexaedra and hyperboloids, and combinations thereof⁶,²⁹. However, all these attempts have in common that they still represent snow as a

<figure><img src="image_5.png" /></figure>

e-mail: alvaro.robledano-perez@univ-grenoble-alpes.fr

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, Grenoble, France. ²Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Etudes de la Neige, Grenoble, France. ³CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France.</footnote>
<figure><img src="image_6.png" /></figure>

**Fig. 1** | From ice crystals to snow microstructure. a Photograph of an ice crystal in the atmosphere (stellar dendrite - K. Libbrecht). b Photograph of natural snow on the ground (decomposing and fragmented precipitation particles - Météo-France). White scale bars: 1 mm. c, d Photographs of natural snow on the ground (rounded grains and depth hoar - F. Dominé). White scale bars: 1 mm. e–g 3D microstructure of three natural snow samples as revealed by X-ray tomography. From left to right:

precipitation particles, rounded grains, and depth hoar. **h, i** Illustration of the absorption enhancement parameter *B* and the geometric asymmetry parameter *g*^G^ in a 2D space. In **h**, the red arrows represent a photon path accounting for refraction and internal reflections, while the white arrows represent propagation in a straight line. In **i**, the white arrows represent an incoming set of photons, and the red arrows represent an ensemble of possible outgoing photon paths.

collection of particles with well-defined shapes. Conversely, Malinka³⁰ considered snow as a two-phase random medium without any particular shape, but it is still unclear whether this approach is fully valid for representing the microstructure of natural snow³¹,³². Despite all these attempts to derive the optical properties of snow by fitting various shapes or hypotheses onto the complexity of the snow microstructure, little is known about its actual role. In simpler words, the “optical shape” of snow—a concept we introduce here following that of the “optical diameter”—remains largely unknown.

Here, we address the question of the optical shape of snow, across the wide diversity of snow microstructures. To do so, we apply a ray-tracing model to simulate the path and attenuation of light in 3D images issued from micro-computed tomography³³⁻³⁵, which over two decades has provided detailed knowledge of the microstructure of snow with a resolution of a few micrometers (Fig. 1e–g). From the simulation results we deduce the values of two main optical shape parameters, namely the absorption enhancement parameter *B* and the geometric asymmetry parameter *g*^G^²⁷ (Fig. 1h, i). *B* quantifies the lengthening of the light path within the absorbing phase (the ice) due to refraction and internal reflections, and is important for accurately estimating the light absorption in snow. *g*^G^ quantifies the tendency of the medium to scatter light forward or backward, and is essential to predict how deep light can penetrate into snow, or conversely how easily it is reflected back to the atmosphere. Both parameters characterize the optical shape independently of the size. Snow optical models, such as those implemented in climate models³⁶⁻³⁸, directly or indirectly rely on prescribed values of these parameters, which are most of the time set to the values for spheres. *B* and *g*^G^ have previously been estimated for individual ice particles with particular shapes²⁷,³⁹⁻⁴² or indirectly estimated from macroscopic measurements on snow⁴³,⁴⁴. However, direct estimation of *B* and *g*^G^ for natural snow along with their variability is still unknown⁴⁵. In the present study, our simulations provide an accurate estimate of the range of these fundamental parameters for natural snow, paving the way for a more realistic representation of snow in optical and climate models.

## Results

### The optical shape of natural snow

The shape parameters *B* and *g*^G^ were computed over 33 snow microstructure images, that cover most of the snow types referenced in the international classification of seasonal snow on the ground⁴ (see Supplementary Table 1). The computations are made with the Rough Surface Ray-Tracing (RSRT) model⁴⁶,⁴⁷, which has been extended to simulate light propagation in 3D microstructure images (see Methods). For every image, simulations are run in the visible and near-infrared (NIR) spectral region, every 50 nm from 400 nm to 1400 nm. The simulations track 10⁶ photons through the microstructure and report the energy reflected back and the profile of energy within the snow-pack, which are then used to deduce *B* and *g*^G^ (macroscopic method hereinafter). In addition, we implemented a microscopic approach (called geometric method) that records the traveled distance within the ice, and the direction changes between entering and exiting the ice (see Methods).

The parameters *B* and *g*^G^ of all snow samples at 900 nm are presented in Fig. 2, along with previous estimates obtained for idealized geometric shapes. Figure 2 also shows the values predicted by the two-phase random medium theory³⁰, where *B* is related to the ice refractive index *n* by $B = n^2$. *g*^G^ can also be expressed in terms of *n*, in a less trivial way (see Supplementary Methods 1). We use here the 900 nm results for a comparison with previous studies⁴⁰,⁴⁸. These results are as well a compromise between the lower wavelengths where the modeling uncertainties are higher (Supplementary Fig. 1) and the higher wavelengths where the assumption of low absorption of ice is less valid than at shorter wavelengths. The absorption enhancement parameter *B* of natural snow clearly clusters around 1.7 (mean of the 33 samples ± 1 standard deviation: 1.68 ± 0.02 in Fig. 2a, 1.70 ± 0.00 in Fig. 2c), while for idealized shapes it spans a larger range between 1.25 for spheres to 1.84 for fractals. However, except fractals, all other featured shapes have a *B* smaller than 1.7, meaning that natural snow absorbs energy more efficiently than these idealized shapes. The value of 1.7 is comparable to that experimentally retrieved in⁴³ (1.6 ± 0.1), and is slightly higher than the one in⁴⁵ (1.49). It matches the predicted value for the
<figure><img src="image_7.png" /><figcaption>**Fig. 2** | The optical shape of snow. a, b Absorption enhancement parameter B and geometric asymmetry parameter $g^G$ (and combinations) of snow at 900 nm, retrieved with the macroscopic method. c, d Idem, retrieved with the geometric method. Note that albedo and light penetration depend on other factors than shape, in particular on grain size, so the representation in b, d must be interpreted at equal snow grain size. In every panel, the dark symbols correspond to geometric</figcaption></figure>

shapes reported in the literature (see Supplementary Table 2) and the two-phase random medium, labeled in a as Malinka (2014) (see Supplementary Methods 1). The colored ones correspond to the 33 natural snow samples, depending on the snow type¹: Precipitation Particles (PP), Decomposing and Fragmented precipitation particles (DF), Faceted Crystals (FC), Depth Hoar (DH), Rounded Grains (RG) and Melt Forms (MF).

two-phase random medium, that is $B = 1.70$ at $\lambda = 900$ nm³. Over the 600 to 1200 nm range, B barely varies, less than 2% with the geometric method and no more than 7% with the macroscopic method (Supplementary Fig. 2). For wavelengths shorter and longer than 600 nm and 1200 nm, respectively, the macroscopic method is less accurate (see Methods) but the geometric method still yields $B = n^2$. Hence, the absorption enhancement of natural snow do not depend on the microstructural details. Instead, B is virtually constant and equal to $n^2$ over the 400 to 1400 nm range.

The geometric asymmetry parameter $g^G$ of natural snow varies between 0.60 and 0.73 (0.65 ± 0.03) at 900 nm with negligible variations across the visible and NIR range (Fig. 2a, c and Supplementary Fig. 2), while for the geometric shapes it spreads from 0.50 for fractals to 0.79 for spheres. This means that snow is less forward-scattering than spheres, but more than most of the featured geometric shapes. These values agree with those presented in previous studies ($g^G$ = 0.68 at $\lambda$ = 900 nm⁴⁵; $g^G$ varying from 0.66 to 0.73⁴⁰ or 0.64 to 0.66 for a more limited dataset⁴⁵).

The albedo and the light penetration depth are controlled by the grain size (via the specific surface area, SSA) and the combination of B and $g^G$. Indeed, the ratio $B/SSA(1-g^G)$ governs the influence of the size and the optical shape on the albedo and the product $B(1-g^G)$SSA on the light penetration⁴⁰. In the representation in Fig. 2b, d, two snowpacks with equal size (i.e. same SSA) but different shapes have the same albedo (respectively light penetration) only if the shapes have the same ordinate (respectively abscissa). In other words, snow samples with equal size but different shape may have different albedo or light penetration. Natural snow spans a region distinct from that of the geometric shapes in the 2D space defined by these quantities (Fig. 2b, d), implying that none of the studied geometric shapes can be used to satisfactorily simulate snow albedo and light penetration at the same time.

Interestingly, a relationship arises between the albedo and the snow type. For a given SSA, fresh snow (PP), like spheres, is a relatively inefficient reflector, while rounded grains (RG) nearly behave as other idealized shapes, such as cylinders (Fig. 2b, d and Supplementary Fig. 3). This is counter-intuitive as rounded grains or melt forms (MF) have the most spherical shape. Regarding light penetration, natural snow behaves similarly to cylinders, and roughly halfway between spheres and fractals (Fig. 2b, d and Supplementary Fig. 3). Fresh snow is, however, more penetrating than the rest of snow types for a given SSA. Even if the spectral albedo of non-spherical shapes can be estimated using spheres by scaling their radius⁴⁹, light penetration depth in a medium with spheres is approximately twice longer than in snow with the same SSA. These results show that, in order to represent these quantities, natural snow should not be represented by the geometric shapes that have been commonly implemented in radiative transfer models, and in particular by spheres.

## Towards a universal representation of snow microstructure in optical models

To further understand why all the featured snow types lead to the same constant value for absorption enhancement, we investigate with the ray-tracing model how the value of B varies when an idealized shape is progressively deformed. To this end we explore three shapes: a sphere, a cube and a convex shape without any symmetry, unlike the first two (Supplementary Methods 2 and Supplementary Fig. 4). For the sphere and the convex shape, the number of triangular facets used to generate the surface is gradually decreased, while for the cube, each of the 8 corners is translated in space in a random direction, with an increasing distance. In all cases, this results in increasingly deformed shapes.

The parameters B and $g^G$ of the three shapes at 900 nm, computed with the geometric method, are presented in Fig. 3. Interestingly, for
<figure><img src="image_8.png" /><figcaption>**Fig. 3 | Optical shape parameters of gradually deformed idealized shapes.**
Variations of the optical shape parameters (absorption enhancement parameter *B* (top); geometric asymmetry parameter *g*<sup>G</sup> (bottom)) of a gradually deformed</figcaption></figure>

sphere (left), cube (center) and analytical convex shape (right). All the simulations
are performed at 900 nm with the geometric method.

the convex shape, *B* is constant, while for spheres and cubes, *B* pro-
gressively converges as the deformation increases to the value found
for natural snow and predicted for the two-phase random medium
(1.7 - Fig. 3, top panel). The theoretical value of *B* for spheres, around
1.25<sup>9,40</sup>, is only obtained for the sphere with the largest number of
facets (≈5 million), suggesting that even the smallest deviation from
this perfection has large consequences for optical properties. This is
relevant to understand that, even if rounded grains or melt forms may
look spherical, their *B* and *g*<sup>G</sup> values considerably differ from those of
"perfect" spheres.

The underlying reason explaining why and when B = n² can be
established from a series of fundamental studies in mathematics,
ecology, optics and nuclear physics⁵⁰⁻⁵³. The absorption within a
weakly-absorbing particle is proportional to the mean path traveled by
photons in the particle, and B measures how this distance is increased
compared to the propagation in a straight line, in the case of diffuse
illumination. B is influenced by two effects, (i) how the photons are
focused as they enter the particle (refraction), and (ii) the mean dis-
tance traveled by photons in the particle. The first effect introduces a
factor n² and is independent of the particle shape as demonstrated
theoretically and experimentally⁵⁴,⁵⁵. The second effect introduces a
factor of exactly 1 (thus leading to B = n²) in several cases: for non-
refractive particles (n = 1) the photons propagate in straight lines and
the mean distance traveled in the particle <|> is given by the Cauchy
formula <|> = 4 V/S (V and S are the volume and surface area of the
particle). The same mean distance is obtained for refractive (n > 1)
particles composed of a scattering material⁵⁵ because the Cauchy
formula holds for a wide class of random walks⁵¹,⁵³. The reason is the
compensating effect of scattering: longer tortuous paths are balanced
by short paths that escape quickly from the particle. However, ice is
not a scattering material. In that case, the mathematical theory of
billiards can be applied to photons bouncing inside a particle⁵⁰, and it
was shown⁵² that if the photons traverse the entire particle in all
directions perfectly uniformly, the mean distance is again given by the
Cauchy formula, which implies B = n². Some billiards (i.e., shapes) are
ergodic and verify this isotropy condition for any refractive index.
Conversely, idealized shapes such as spheres and cubes, are non-
ergodic and some regions may not be uniformly explored by photons
coming from the outside, especially if the refractive index is larger than
a shape-dependent critical value⁵². As these unexplored regions

generally correspond to very long paths that are only accessible
through internal scattering, the mean traveled distance decreases,
leading to *B* < *n*² as observed in Fig. 4 for spheres and cubes. Note also
that strong absorption also reduces the very long paths, leading to a
decreased *B* (Supplementary Fig. 5). To conclude, the fact that we find
*B* = *n*² for all the investigated snow samples in the visible and NIR
spectral region strongly suggests that the snow microstructure is
fundamentally ergodic.

To investigate whether this result applies to materials other than
snow, we computed *B* for different *n* values with the ray-tracing model.
We find that the value of *B* for near-perfect spheres closely follows the
analytical expression for spheres⁹, while the equality between *B* and *n*²
stands for less symmetric spheres (Fig. 4a), as well as for the large
diversity of snow samples studied here (Fig. 4b). This equality is more
general and is actually an experimental evidence that, in terms of
absorption enhancement, a weakly absorbing porous material like
snow can be represented as a collection of convex particles without
symmetry or as a two-phase random medium. Representing such
porous media in these ways may be useful and crucial for refining the
computation of other shape-dependent optical properties⁵⁶, with a
wide range of applications well beyond the snow optics community,
such as the optical characterization of pharmaceutical powders⁵⁷ or
solar cell design⁵⁸.

In contrast to *B*, *g*<sup>*G*</sup> is almost unaffected by the shape deformation (Fig. 3, bottom panel), although it is more sensitive to the particle shape or the type of snow (Fig. 2). For spheres, *g*<sup>*G*</sup> closely matches the value obtained by theoretical calculations<sup>9</sup> and in other studies<sup>27,40</sup>. For cubes, *g*<sup>*G*</sup> decreases very slightly when the random displacement is applied to the corners, breaking the symmetries between the faces. While *B* is very sensitive to the rare and long paths discussed above, we conclude that *g*<sup>*G*</sup> is more sensitive to the few internal reflections experienced by the photons, which determine the ability of snow to scatter light forward. Consequently, *g*<sup>*G*</sup> does not have a unique value for snow, but mainly spans the interval 0.62–0.68, with a slightly higher value for fresh snow. This range is however much smaller than that estimated from idealized shapes, and importantly, it does not contain the value for spheres (0.79), still commonly implemented in some climate models<sup>36,37</sup>. The values of the two-phase random medium and the convex particle (0.67 and 0.63, respectively) are more consistent but they do not represent the full range of values across the
<figure><img src="image_9.png" /></figure>

**Fig. 4** | Variations of the absorption enhancement parameter *B* with a varying ice refractive index *n*. *a* Variations for spheres. The number of triangular facets used to generate the surface is gradually increased. The BB 1974 curve is the analytical expression for spheres of *B* in terms of *n* (Eq. 9 in¹). *b* Variations for natural snow. 9 samples are considered, covering most of the main snow types. Three 3D

images of snow samples are shown to illustrate the wide diversity of snow micro-
structures. The simulations are done with the geometric method, and the ice
absorption coefficient γ is kept constant and corresponding to the ice absorption
at 900 nm.

diverse types of snow. Even if different *B* and *g*<sup>*G*</sup> values might be
expected for under-represented snow types in our study, such as very
peculiar samples of fresh snow or strongly developed depth hoar, our
results show that the behavior of snow in the 400 to 1400 nm range is
less variable and less shape-dependent than thought so far. This con-
stitutes a step towards a more accurate and universal representation of
snow in optical models.

Discussion

Climate simulations suggest that an ambiguous treatment of snow
shape leads to large uncertainties in the estimation of surface radiative
forcing⁶,⁷ and global air temperature of up to 1.17 K⁸. This uncertainty is
due to the high sensitivity of the global temperature to the snow
albedo, enhanced by potent climate feedbacks²⁹⁻⁶¹. These uncertainties
in snow albedo are driven by *B* and *g*<sup>*G*</sup> through the single-scattering
properties of snow. While these parameters might not be directly used
in such relevant models, it is possible to convert these quantities into
other optical properties, such as the asymmetry parameter *g* and the
single-scattering albedo *ω*. *g* determines the first order angular varia-
tions of the phase function and is directly related to the previously
studied *g*<sup>*G*</sup>. The single-scattering co-albedo (1−*ω*) depends linearly, for
a given particle size, on *B* (see Methods). Only recently snow radiative
transfer schemes implemented in climate models, such as SNICAR-
ADv3, started to consider non-spherical shapes for *g*, with para-
meterizations based on spheroids, hexagonal plates and fractals<sup>6,7,38,62</sup>.
Still, the spherical assumption is used for *B* (indirectly via *ω*). This
simplification has a direct impact on the albedo, and also strong con-
sequences for the light penetration depth, as spheres are roughly twice
more penetrating than natural snow (Fig. 2). This has implications for
the thermal regime of the snowpack<sup>63</sup> and for the transmittance of
snow over ground<sup>45</sup> or over sea ice<sup>26</sup>. To our knowledge, only in<sup>s</sup> non-
spherical values of *B* and *g*<sup>*G*</sup> are considered. Their 1.17 K global air
temperature change was obtained by varying together *g* from 0.89 to
0.78 and *B* from 1.25 to ≈ 1.62 (indirectly via *ω*<sup>64</sup>).

Here we provide strong observational constraints on *B* and *g* for
snow, which shall help reduce uncertainties in climate studies by
successfully simulating snow albedo and light penetration at the same

time. *B* is universal and equal to *n*², and according to our simulations,
the value of *g* for snow over the 33 samples is *g* = 0.82 ± 0.01, with
negligible variations in the 400 to 1400 nm wavelength range.
Although performing climate simulations with these updated values is
beyond the scope of this study, we propose a simple estimation of the
impact on global air temperature, based on the study from⁸ and by
considering the quantities represented in Fig. 2b, d (see Methods). By
using the constant values of *B* = 1.7 and *g* = 0.82 (or equivalently
*g*<sup>*G*</sup> = 0.65) instead of the values for the non-spherical shape, the simu-
lated global annual-mean air temperature would shift by roughly 0.6 K,
assuming similar sensitivity of temperature to snow albedo as in⁸.
Moreover, the narrow range of values found here for natural snow
would drastically reduce the uncertainties due to the equivocal impact
of snow morphology, dropping from 1.17 K to approximately 0.4 K.
Beyond climate simulations, the refined knowledge of the optical
shape of snow obtained in this study will be beneficial wherever snow
optics matters, from snow photochemistry to remote sensing algo-
rithms, solving the long-standing issue of the optical shape of snow.

Methods

Snow microstructure from X-ray tomography

The 3D snow microstructure images have been acquired over the last
decades by the Centre d'Études de la Neige (Météo-France - CNRS). The
dataset presented here (33 images) includes mainly seasonal snow
collected in the French Alps, as well as several samples resulting from a
series of laboratory metamorphism experiments³²,⁶⁵⁻⁶⁹. The X-ray
tomography was performed at various resolutions, from 4.91 µm to
11.65 µm. Depending on the experimental set-up, the sample volume
varied between 2.5 × 2.5 × 2.5 and 9.5 × 9.5 × 9.5 mm³. For image
acquisition, samples were first impregnated with liquid
1-chloronaphthalene around -10 °C, and placed below -20 °C, thus
forcing the chloronaphthalene to freeze. Small cylinders were then
machined at -30 °C, inserted in plexiglas caps and fixed on copper
columns for their later insertion into CellStat, a cold cell specifically
designed for tomography of frozen samples at ambient
temperature⁶⁸,⁷⁰. After tomography, the differences in X-ray attenua-
tion between the ice, 1-chloronaphthalene and remaining air bubbles
were used to segment the reconstructed 3D grey-level images into 3D binary images. Surface meshes were then obtained through an automated pipeline using iso2mesh<sup>71</sup>. For each image, the function 'v2s' of this mesh generation toolbox was used with the option 'cgalmesh', ensuring the most robust path to product meshes from binary volumes. The maximal deviation from the 0.5 isosurface was set below 0.5 and the maximum radius of the Delaunay spheres was chosen smaller than 15 pixels. For each closed surface, the mesh was then automatically checked for orientation consistency of the facets. A final manual verification and correction step was realized using MeshLab<sup>72</sup>.

The whole description of the dataset (type of snow, resolution, etc) is available in the Supplementary Table 1. The corresponding snow microstructural properties (density, SSA) were computed directly over the mesh with the trimesh Python package<sup>73</sup>. Visualizations of several of these images are available in Fig. 1e–g.

## The ray-tracing model

The existing Rough Surface Ray-Tracing (RSRT) model, originally designed to simulate albedo over rough surfaces and topography<sup>46,47</sup>, has been adapted to trace light propagation in 3D microstructure images. In this Monte Carlo photon-tracking approach, a set of photons (rays) is followed through the snow from their source to termination (absorption or escape at the top of the sample), relying on the geometric optics approximation (i.e., we consider the microstructure features to be much larger than the wavelength) to simulate the ray path. The light propagation is governed by both absorption within the ice phase (which is wavelength-dependent) and the fundamental laws of reflection and refraction at each air–ice interface, as in<sup>74</sup> and<sup>16</sup>. The origin and initial direction of each ray is randomly generated above the snow microstructure image. When a ray encounters an ice–air interface, the choice between reflection and refraction is random, and depends on the Fresnel’s law of reflectance (Supplementary Methods 3). The chosen ray carries all the incident energy. When traveling through the ice phase, part of the energy carried by the ray is lost by absorption, following an exponential decay proportional to the traveled distance in the ice phase (i.e. Beer’s law). This is wavelength-dependent as it considers the ice absorption coefficient γ:

$$y = \frac{4\pi}{\lambda} n_i(\lambda) \qquad (1)$$

where $n_i(\lambda)$ is the imaginary part of the ice refractive index<sup>75</sup>.

To simulate a semi-infinite snowpack, a replication algorithm is applied when a ray goes through the boundaries of the original snow microstructure image. Replicating periodically the source image allows simulating a macroscopic snowpack from a single snow microstructure image, largely reducing the computational cost of this approach (in particular with respect to the memory storage of the mesh). It also has some drawbacks, as we consider the snowpack to be homogeneous and single-layer, with snow properties defined by the source microstructure image, such as density or specific surface area (SSA). However, this does not prevent the use of the intrinsic B and g<sup>G</sup> quantities in a multi-layered snowpack, where the properties vary vertically<sup>76</sup>.

Every single ray is traced until it escapes the simulated semi-infinite snowpack in the upward direction, or until its energy goes below a defined, very low threshold (we consider the ray to be absorbed). N = 10<sup>6</sup> rays are launched in each simulation to reach a reasonable accuracy (uncertainty in this Monte-Carlo framework decreases as 1/√N).

## B and g<sup>G</sup> computation: the macroscopic method

Two different methods have been implemented to derive the optical shape parameters, B and g<sup>G</sup>, from 3D images. The hereinafter called

macroscopic method relies on the reflected energy by the snowpack and the vertical profile of energy in the snowpack. Starting from the approximate asymptotic radiative transfer (AART) theory<sup>27</sup>, and with the formalism developed in<sup>40</sup>, the bi-hemispherical albedo α (hereafter albedo) and the asymptotic flux extinction coefficient $k_e$ are expressed in terms of the optical shape parameters B and g<sup>G</sup> by:

$$\alpha \approx \exp \left( -4 \sqrt{\frac{B\gamma V}{3\Sigma(1-g^G)}} \right) \qquad (2)$$

$$k_e \approx \frac{\rho}{\rho_{\text{ice}}} \sqrt{\frac{3B\gamma\Sigma}{V}(1-g^G)} \qquad (3)$$

from where we formulate B and g<sup>G</sup> as:

$$B \approx \frac{-\rho_{\text{ice}}k_e(\lambda)\ln(\alpha(\lambda))}{4\rho\gamma(\lambda)} \qquad (4)$$

$$g^G \approx 1 + \frac{16k_e(\lambda)}{3\rho\text{SSA}\ln(\alpha(\lambda))} \qquad (5)$$

using the snow SSA = 4Σ / (Vρ<sub>ice</sub>), ρ the snow density and ρ<sub>ice</sub> the ice density (i.e. 917 kg m<sup>-3</sup>). Σ and V are, respectively, the average projected area and the average volume of a particle. Albedo differences of 0.002 are found between a Mie+DISORT model (e.g.<sup>77</sup>) – a robust and accurate radiative transfer model for spherical particles – and the AART theory using B and g<sup>G</sup> relevant for spheres (1.25 and 0.79, respectively), and a SSA of 16.4 m<sup>2</sup> kg<sup>-1</sup>. These albedo differences yield a relative error of ≈ 1% in the (B, g<sup>G</sup>) computations (≈ 0.02 and ≈ 0.007, respectively), which we consider negligible. These expressions remain thus valid in the limit of low absorption, which is globally true in the visible and NIR spectral region (400 to 1400 nm).

In the RSRT model we compute the albedo α as the ratio between the number of photons escaping the semi-infinite snowpack with an upward direction and the total number of photons launched. Using a collection of monodisperse spherical particles, the spectral relative error between albedos computed with RSRT and those predicted by the Mie+DISORT model was lower than 3% in the 400 to 1200 nm wavelength range (slightly higher in the 1200 to 1400 nm range - Supplementary Fig. 6). The asymptotic flux extinction coefficient $k_e$ is computed counting the intensity carried by the photons passing through a given horizontal plane z. This intensity shows an exponential decrease with depth, from where we fit a Beer-Lambert law (i.e. I(z) ≈ I(z = 0) exp(-k<sub>e</sub> z)) to obtain $k_e$ and eventually compute B and g<sup>G</sup> with Eqs. 4, 5. The modeling uncertainty of α and $k_e$ is treated here with a Bayesian framework.

## Bayesian treatment of uncertainties

To ensure the accuracy of the computed optical shape parameters with the macroscopic method, we implemented a Bayesian model to quantify the uncertainties. Two optical quantities are derived from the ray-tracing simulations: (i) the albedo, and (ii) the flux extinction coefficient $k_e$. These unknown quantities are then described using some known probability distributions (named priors), which are updated using Bayes’ theorem, a process called inference. The resulting posterior distribution provides an estimation of the albedo and $k_e$ modeling uncertainties<sup>78</sup>.

Here, for each simulation we describe the prior distribution of albedo with a normal distribution N (μ, σ<sup>2</sup>), with μ being the computed albedo and σ = 1/√N. For $k_e$, it is less direct. The probability of finding a photon exponentially decreases with depth z, so we could describe this as a Bernoulli process. Consequently, the prior distribution of
observing *n* photons at a certain depth *z* is given by a Binomial distribution $B(n, p)$, with $p = \exp(-k_e z)$.

The inference is then performed by means of a Python package⁷⁹. To obtain the posterior estimates, the model fitting is based on samples drawn from the posterior distribution using Markov Chain Monte Carlo (MCMC) methods. In particular, a No-U-Turn Sampler (NUTS) is implemented here (8000 samples). Instead of using only the maximum likelihood estimation, all the posterior estimations of both albedo and $k_e$ are then introduced in Eqs. 4, 5 to obtain a full distribution of the optical shape parameters *B* and $g^G$ (Supplementary Fig. 1).

## B and $g^G$ computation: the geometric method

We also implemented a more direct approach (called geometric method) where the geometric definition of the shape parameters is computed by means of ray-tracing from the actual snow microstructure. This method relies only on the geometric optics approximation. The absorption enhancement parameter *B* can be defined by the lengthening of the photon path in the ice phase due to internal reflections and refraction with respect to strictly straight lines⁵⁶. The asymmetry parameter (i.e. $g^G$ in this geometrical optics framework) is defined by the scattering phase function⁸⁰,⁸¹, and is expressed as:

$$g^G = \frac{1}{2} \int_{-1}^{1} d(\cos \Theta) \cos \Theta p(\Theta) \quad (6)$$

with $\Theta$ the scattering angle and $p(\Theta)$ the scattering phase function, normalized as:

$$\frac{1}{2} \int_{-1}^{1} d(\cos \Theta) p(\Theta) = 1 \quad (7)$$

In the RSRT model, in order to compute the scattering phase function, we record the incident and outbound ray direction when entering and going out the ice phase, respectively. The geometric asymmetry parameter $g^G$ is directly deduced from Eq. 6.

## Relating the optical shape parameters to the snow single-scattering properties

The two optical shape parameters, *B* and $g^G$, can be used in climate modeling. The first step would be to use them to calculate other fundamental snow optical properties, such as the single-scattering albedo $\omega$ and the asymmetry parameter *g*. Defining the single-scattering albedo (1 - $\omega$) as the ratio of absorption to extinction coefficients, and following the formalism developed in⁴⁰ (Eqs. (1–6)), we can relate *B* to $\omega$ by:

$$(1 - \omega) = B \gamma \frac{V}{2\Sigma} \quad (8)$$

where the snow SSA can be introduced by using SSA = $4\Sigma / (V\rho_{ice})$:

$$(1 - \omega) = B \gamma \frac{2}{SSA \rho_{ice}} \quad (9)$$

The asymmetry parameter *g* is simply the average of the geometric and the diffraction terms ($g^G$ and $g^D$, respectively). In this framework, where the wavelength is small enough compared to the particles, diffraction is mainly forward (i.e. $g^D \approx 1$), so that:

$$g = \frac{g^G + 1}{2} \quad (10)$$

For more advanced models, *g* can be used to parameterize the phase function. For instance in DISORT⁸², that uses the Legendre

polynomial decomposition of the phase function, *g* appears to be the coefficient of the first order polynomial.

For models that only require the broadband value of the snow albedo, the values of *B* and $g^G$ found here for natural snow could be used in snow radiative transfer models that rely on such parameters (e.g. TARTES⁴⁰) to eventually derive an updated snow albedo parameterization⁸³,⁸⁴.

## Estimation of temperature uncertainty reduction

In⁸, it was found that if the shape spans the range from spheres to the Optimized Habit Combination (OHC), the global annual-mean air temperature varies by 1.17 K. Since the snow albedo depends on the ratio $\Gamma = B / (1 - g^G)$, in order to estimate the reduction of the uncertainties related to the optical shape of snow in climate modeling, this quantity and its variations for natural snow are evaluated and compared to those in⁸.

Using the values for the sphere ($B = 1.25$, $g^G = 0.79$) and the values for the OHC ($B \approx 1.62$, $g^G \approx 0.56$)⁶⁴, the $\Gamma$ range explored by⁸ is:

$$\Gamma_{\text{Räisänen}} = \frac{B_{\text{sph}}}{1 - g_{\text{sph}}^G} - \frac{B_{\text{OHC}}}{1 - g_{\text{OHC}}^G} = 2.27 \quad (11)$$

In the present study, we obtained a reduced range for natural snow. *B* can be considered constant and equal to 1.7, and $g^G$ varies mainly between 0.62 and 0.68 (0.65 ± 0.03 - Fig. 2). If we use our values instead of those for the sphere and the OHC, the $\Gamma$ range explored in this study is:

$$\Gamma_{RSRT} = \frac{B_{\text{snow}}}{1 - g_{\text{snow,upper}}^G} - \frac{B_{\text{snow}}}{1 - g_{\text{snow,lower}}^G} = 0.84 \quad (12)$$

From these values we conclude that the shape uncertainty is reduced by a factor $\Gamma_{Räisänen}/\Gamma_{RSRT} \approx 3$. Assuming linear sensitivity of global temperature to albedo and linear dependency between albedo and $\Gamma$ (valid for small perturbations), the temperature uncertainty is reduced then from 1.17 K to 0.43 K.

## Method limitations

### Snow microstructure.

(1) Our results are based on a finite set of snow images. Although meant to be representative of the diversity of snow, it is limited to 33 images, so that different *B* and *g* values cannot be ruled out for very peculiar, under-represented snow morphologies.

(2) Note that only pure snow, i.e. without light-absorbing particles (LAPs), is considered in this study. This is due to the fact that the size distribution of dust particles is much lower than the resolution limit of our X-ray tomography images (a few µm). A recent study showed the feasibility of capturing the motion of dust particle aggregates in dry snow⁸⁵, which could potentially open the way to determine the impact of LAPs on the optical properties of snow by ray-tracing. Work is underway to implement this feature in our ray-tracing model.

### Ray-tracing model.

(1) Wave properties of light (diffraction and polarization) are not considered in our ray-tracing model. However, several studies showed that these geometric optics simplifications do not prevent to correct simulate snow reflectance¹⁶,³⁰,⁷⁴,⁸⁶.

(2) Some numerical errors may occur with the most complex snow microstructures (mainly corresponding to fresh snow). This is due to the resolution limit of our current imaging technique, that fails to resolve the most detailed features of such microstructures. This leads to slight artifacts in the mesh generation process and in turn may induce some errors in the photon-tracking method. This is particularly true for the longest photon
paths, at the shortest wavelengths, when the ice absorption is extremely weak.

(3) The definition of a single-scattering event in this study slightly differs from what is typically assumed for unconnected particles⁴¹. Here, a single-scattering event ends when the ray first exits the ice phase (including reflection at the entrance), whereas in the common definition a ray may enter and exit several times if the particle is concave, before finally escaping the particle. Unfortunately, this common definition requires extracting independent particles from 3D images, which is somewhat arbitrary, since the ice phase is usually mostly connected. Individual snow grains can however be defined as zones separated by regions of potential mechanical weakness (e.g.⁸⁷) though these individual snow grains are still connected by ice. The surface area of these ice-ice contacts is nevertheless small compared to the ice-air interface area⁶⁶. In conclusion, for natural snow, and in order to use the snow microstructure images as is, this uncommon but pragmatic definition of a single-scattering event was used here in the geometric method to derive $g^G$.

### Model uncertainties.

(1) The $B$ and $g^G$ values presented in the text and figures with the macroscopic method correspond to the mean values of the resulting $B$ and $g^G$ distributions. Their dispersion is usually small, with a mean standard deviation over the whole dataset and the 600–1200 nm wavelength range $\bar{\sigma}$ of 0.02 and 0.004 for $B$ and $g^G$ estimations, respectively. For $B$ it is therefore unlikely that the conclusion $B=n^2$ is affected (Supplementary Fig. 1). The same applies to the estimation uncertainty of $g^G$, which is also considerably smaller than the differences between natural snow and the considered geometric shapes.

(2) In relation with the precedent limitation about the modeling uncertainties, the macroscopic method is less accurate at the shortest and longest wavelength range (400–600 nm and 1200–1400 nm, respectively). Below 600 nm, the albedo computation needs to be extremely accurate to derive a precise estimate of the ($B$, $g^G$) parameters, which is very computation-costly with a Monte-Carlo approach. This happens because the snow albedo $\alpha$ in this spectral region is close to 1, and in the limit of $\alpha \approx 1$, the In ($\alpha$) dependence becomes very close to zero and in particular, a slight underestimation of the albedo (likely due to the numerical cutoff of the photons) leads to a large overestimation of $B$ (Supplementary Fig. 1). Above 1200 nm, the limitation comes from the underlying asymptotic radiative transfer theory that is only valid in the low absorption limit, which might not be fully respected at these longer wavelengths. However, it is important to note that the geometric method still yields $B=n^2$ below 600 nm and above 1200 nm. For longer wavelengths, the geometric optics approximation (where particles need to have dimensions much larger than the incident wavelength), limits as well the geometric method. To summarize, our results are valid over the 400 to 1400 nm range, which encompasses the most important part of the solar spectrum ($\approx$ 85% of the solar irradiance at the surface is within this range).

(3) In the macroscopic method, some very small uncertainties in the ($B$, $g^G$) computation might come from the input macroscopic quantities (density and SSA). In particular, $B$ relies on the snow density $\rho$ (Methods, Eq. 4), and $g^G$ relies on both $\rho$ and SSA (Methods, Eq. 5). These quantities can be estimated with a 2% accuracy⁸⁸. Propagating this uncertainty into the equations, this is equivalent to errors of 0.03 and 0.02 in $B$ and $g^G$, respectively. Moreover, our computed $\rho$ and SSA values compare well to measurements over the voxelized microstructure images in³², with similar accuracy. However, it is important to note that the geometric method, which does not suffer from these

uncertainties in the input macroscopic variables, yields equivalent results to the macroscopic method.

### Data availability

The generated geometric shapes have been deposited in the PerSCIDO platform and are available from https://doi.org/10.18709/perscido.2023.06 ds392. Source data are provided in this paper.

### Code availability

The simulation results and codes to generate the figures have been deposited in the PerSCIDO platform and are available from https://doi.org/10.18709/perscido.2023.06 ds392. The ice refractive index is computed with the snowoptics library available from https://github.com/ghislainp/snowoptics (last access: 16 June 2023) and https://doi.org/10.5281/zenodo.3742138.

### References

1. Lawson, R. P. et al. Deployment of a tethered-balloon system for microphysics and radiative measurements in mixed-phase clouds at ny-ålesund and south pole. *J. Atmos. Oceanic Technol.* **28**, 656–670 (2011).

2. Greenler, R. *Rainbows, Halos, and Glories* (Cambridge U. Press, Cambridge, UK, 1980).

3. Moilanen, J. & Gritsevich, M. Light scattering by airborne ice crystals – an inventory of atmospheric halos. *J. Quant. Spectrosc. Radiat. Transfer* **290**, https://doi.org/10.1016/J.JQSRT.2022.108313 (2022).

4. Fierz, C. et al. The international classification for seasonal snow on the ground (UNESCO/IHP, Paris, France, 2009).

5. Pinzer, B. R., Schneebeli, M. & Kaempfer, T. U. Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography. *Cryosphere* **6**, 1141–1155 (2012).

6. He, C. et al. Black carbon-induced snow albedo reduction over the tibetan plateau: uncertainties from snow grain shape and aerosol-snow mixing state based on an updated snicar model. *Atmos. Chem. Phys.* **18**, 11507–11527 (2018).

7. Hao, D. et al. Improving snow albedo modeling in the e3sm land model (version 2.0) and assessing its impacts on snow and surface fluxes over the tibetan plateau. *Geosci. Model Dev.* **16**, 75–94 (2023).

8. Räisänen, P., Makkonen, R., Kirkevåg, A. & Debernard, J. B. Effects of snow grain shape on climate simulations: sensitivity tests with the norwegian earth system model. *Cryosphere* **11**, 2919–2942 (2017).

9. Bohren, C. F. & Barkstrom, B. R. Theory of the optical properties of snow. *J. Geophys. Res.* (1896–1977) **79**, 4527–4535 (1974).

10. Warren, S. G. Optical properties of snow. *Rev. Geophys. Space Phys.* **20**, 67–89 (1982).

11. Liou, K. N. & Takano, Y. Light scattering by nonspherical particles: remote sensing and climatic implications. *Atmos. Res.* **31**, 271–298 (1994).

12. Grenfell, T. C. & Warren, S. G. Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation. *J. Geophys. Res.* **104**, 31697–31710 (1999).

13. Domine, F. et al. Correlation between the specific surface area and the short wave infrared (SWIR) reflectance of snow. *Cold Reg. Sci. Technol.* **46**, 60–68 (2006).

14. Neshyba, S. P., Grenfell, T. C. & Warren, S. G. Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation: 2. hexagonal columns and plates. *J. Geophys. Res.* **108**, 4448–+ (2003).

15. Grenfell, T. C., Neshyba, S. P. & Warren, S. G. Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation: 3. hollow columns and plates. *J. Geophys. Res.* **110**, 17203–+ (2005).
16. Picard, G., Arnaud, L., Domine, F. & Fily, M. Determining snow specific surface area from near-infrared reflectance measurements: numerical study of the influence of grain shape. *Cold Reg. Sci. Technol.* **56**, 10–17 (2009).

17. Stanton, B., Miller, D., Adams, E. & Shaw, J. A. Bidirectional-reflectance measurements for various snow crystal morphologies. *Cold Reg. Sci. Technol.* **124**, 110–117 (2016).

18. Dumont, M. et al. High-accuracy measurements of snow bidirectional reflectance distribution function at visible and NIR wavelengths - comparison with modelling results. *Atmos. Chem. Phys.* **10**, 2507–2520 (2010).

19. Sergent, C., Chevrand, P., Lafeuille, J. & Marbouty, D. Caractérisation optique de différents types de neige. extinction de la lumière dans la neige. *J. Phys. Colloques* **48**, 361–367 (1987).

20. Meirold-Mautner, I. & Lehning, M. Measurements and model calculations of the solar shortwave fluxes in snow on summit, greenland. *Ann. Glaciol.* **38**, 279–284 (2004).

21. Leroux, C., Lenoble, J., Brogniez, G., Hovenier, J. W. & De Haan, J. F. A model for the bidirectional polarized reflectance of snow. *J. Quant. Spectrosc. Radiat. Transfer* **61**, 273–285 (1999).

22. Schaepman-Strub, G., Schaepman, M. E., Painter, T. H., Dangel, S. & Martonchik, J. V. Reflectance quantities in optical remote sensing—definitions and case studies. *Remote Sens. Environ.* **103**, 27–42 (2006).

23. Domine, F. et al. Snow physics as relevant to snow photochemistry. *Atmos. Chem. Phys.* **8**, 171–208 (2008).

24. Erbland, J. et al. Air-snow transfer of nitrate on the east antarctic plateau – part 2: An isotopic model for the interpretation of deep ice-core records. *Atmos. Chem. Phys.* **15**, 12079–12113 (2015).

25. Winton, V. H. L. et al. Deposition, recycling, and archival of nitrate stable isotopes between the air-snow interface: comparison between dronning maud land and dome c, antarctica. *Atmos. Chem. Phys.* **20**, 5861–5885 (2020).

26. Vérin, G., Domine, F., Babin, M., Picard, G. & Arnaud, L. Metamorphism of snow on arctic sea ice during the melt season: impact on spectral albedo and radiative fluxes through snow. *Cryosphere* **16**, 3431–3449 (2022).

27. Kokhanovsky, A. A. & Zege, E. P. Scattering optics of snow. *Appl. Opt.* **43**, 1589–1602 (2004).

28. Ishimoto, H. et al. Snow particles extracted from x-ray computed microtomography imagery and their single-scattering properties. *J. Quant. Spectrosc. Radiat. Transfer* **209**, 113–128 (2018).

29. Saito, M., Yang, P., Loeb, N. G. & Kato, S. A novel parameterization of snow albedo based on a two-layer snow model with a mixture of grain habits. *J. Atmos. Sci* **76**, 1419–1436 (2019).

30. Malinka, A. V. Light scattering in porous materials: geometrical optics and stereological approach. *J. Quant. Spectrosc. Radiat. Transfer* **141**, 14–23 (2014).

31. Malinka, A., Zege, E., Heygster, G. & Istomina, L. Reflective properties of white sea ice and snow. *Cryosphere* **10**, 2541–2557 (2016).

32. Dumont, M. et al. Experimental and model-based investigation of the links between snow bidirectional reflectance and snow microstructure. *Cryosphere* **15**, 3921–3948 (2021).

33. Brzoska, J. B. et al. 3D visualization of snow samples by microtomography at low temperature. *ESRF Newsletter* **32**, 22–23 (1999).

34. Coléou, C., Lesaffre, B., Brzoska, J. B., Ludwig, W. & Boller, E. Three-dimensional snow images by X-ray microtomography. *Ann. Glaciol.* **32**, 75–81 (2001).

35. Schneebeli, M. & Sokratov, S. A. Tomography of temperature gradient metamorphism of snow and associated changes in heat conductivity. *Hydrol. Processes* **18**, 3655–3665 (2004).

36. Holland, M. M., Bailey, D. A., Briegleb, B. P., Light, B. & Hunke, E. Improved sea ice shortwave radiation physics in ccsm4: the impact

of melt ponds and aerosols on arctic sea ice. *J. Clim.* **25**, 1413–1430 (2012).

37. Flanner, M. G., Zender, C. S., Randerson, J. T. & Rasch, P. J. Present-day climate forcing and response from black carbon in snow. *J. Geophys. Res.* **112**, https://doi.org/10.1029/2006jd008003 (2007).

38. Flanner, M. G. et al. Snicar-adv3: a community tool for modeling spectral snow albedo. *Geosci. Model Dev.* **14**, 7673–7704 (2021).

39. Kokhanovsky, A. A. & Macke, A. Integral light-scattering and absorption characteristics of large, nonspherical particles. *Appl. Opt.* **36**, 8785–8790 (1997).

40. Libois, Q. et al. Influence of grain shape on light penetration in snow. *Cryosphere* **7**, 1803–1818 (2013).

41. Saito, M. & Yang, P. Oriented ice crystals: a single-scattering property database for applications to lidar and optical phenomenon simulations. *J. Atmos. Sci* **76**, 2635–2652 (2019).

42. Xu, G., Schnaiter, M. & Järvinen, E. Accurate retrieval of asymmetry parameter for large and complex ice crystals from in-situ polar nephelometer measurements. *J. Geophys. Res.: Atmos.* **127**, https://doi.org/10.1029/2021JD036071 (2022).

43. Libois, Q. et al. Experimental determination of the absorption enhancement parameter of snow. *J. Glaciol.* **60**, 714–724 (2014).

44. Ottaviani, M., van Diedenhoven, B. & Cairns, B. Photopolarimetric retrievals of snow properties. *Cryosphere* **9**, 1933–1942 (2015).

45. Letcher, T., Parno, J., Courville, Z., Farnsworth, L. & Olivier, J. A generalized photon-tracking approach to simulate spectral snow albedo and transmittance using X-ray microtomography and geometric optics. *Cryosphere* **16**, 4343–4361 (2022).

46. Larue, F. et al. Snow albedo sensitivity to macroscopic surface roughness using a new ray-tracing model. *Cryosphere* **14**, 1651–1672 (2020).

47. Robledano, A., Picard, G., Arnaud, L., Larue, F. & Ollivier, I. Modeling surface temperature and radiation budget of snow-covered complex terrain. *Cryosphere* **16**, 559–579 (2022).

48. Krol, Q. & Löwe, H. Relating optical and microwave grain metrics of snow: the relevance of grain shape. *Cryosphere* **10**, 2847–2863 (2016).

49. Dang, C., Fu, Q. & Warren, S. G. Effect of snow grain shape on snow albedo. *J. Atmos. Sci.* **73**, 3573–3583 (2016).

50. Szász, D. Multidimensional hyperbolic billiards. *1701.02955* (2017).

51. Blanco, S. & Fournier, R. An invariance property of diffusive random walks. *Europhys. Lett.* **61**, 168 (2003).

52. Majic, M., Somerville, W. R. C. & Le Ru, E. C. Mean path length inside non-scattering refractive objects. *Phys. Rev. A* **103**, L031502 (2021).

53. Reuss, P. Cauchy's theorem and generalization. *EPJ Nucl. Sci. Technol.* **4**, 50 (2018).

54. Yablonovitch, E. Statistical ray optics. *J. Opt. Soc. Am.* **72**, 899–907 (1982).

55. Savo, R. et al. Observation of mean path length invariance in light-scattering media. *Science* **358**, 765–768 (2017).

56. Libois, Q., Lévesque-Desrosiers, F., Lambert-Girard, S., Thibault, S. & Domine, F. Optical porosimetry of weakly absorbing porous materials. *Opt. Express* **27**, 22983–22993 (2019).

57. Scheibelhofer, O., Wahl, P. R., Larchevêque, B., Chauchard, F. & Khinast, J. G. Spatially resolved spectral powder analysis: experiments and modeling. *Appl. Spectrosc.* **72**, 521–534 (2018).

58. Sychugov, I. Geometry effects on luminescence solar concentrator efficiency: analytical treatment. *Appl. Opt.* **59**, 5715–5722 (2020).

59. Qu, X. & Hall, A. What controls the strength of snow-albedo feedback. *J. Clim* **20**, 3971–3981 (2007).

60. Flanner, M. G., Shell, K. M., Barlage, M., Perovich, D. K. & Tschudi, M. A. Radiative forcing and albedo feedback from the northern
hemisphere cryosphere between 1979 and 2008. Nat. Geosci. **4**, 151–155 (2011).

61. Riihelä, A., Bright, R. M. & Anttila, K. Recent strengthening of snow and ice albedo feedback driven by antarctic sea-ice loss. Nat. Geosci. **14**, 832–836 (2021).

62. He, C. et al. Impact of snow grain shape and black carbon–snow internal mixing on snow optical properties: parameterizations for climate models. J. Clim. **30**, 10019–10036 (2017).

63. Kuipers Munneke, P. et al. The role of radiation penetration in the energy budget of the snowpack at summit, greenland. Cryosphere **3**, 155–165 (2009).

64. Räisänen, P., Kokhanovsky, A., Guyot, G., Jourdan, O. & Nousiainen, T. Parameterization of single-scattering properties of snow. Cryo-sphere **9**, 1277–1301 (2015).

65. Flin, F., Brzoska, J. B., Lesaffre, B., Coléou, C. & Pieritz, R. A. Three- dimensional geometric measurements of snow microstructural evolution under isothermal conditions. Ann. Glaciol. **38**, 39–44 (2004).

66. Flin, F. et al. Furukawa, Y. (ed.) On the computations of specific surface area and specific grain contact area from snow 3D images. (ed. Furukawa, Y.) Proceedings of the 12th International Conference on the Physics and Chemistry (PCI 2010) of Ice held at Sapporo, Japan, on 5–10 September 2010, 321–328 (Hokkaido University Press, Sapporo, Japan, 2011).

67. Calonne, N. et al. 3-D image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy. Cryosphere **6**, 939–951 (2012).

68. Calonne, N., Flin, F., Geindreau, C., Lesaffre, B. & Rolland du Ros-coat, S. Study of a temperature gradient metamorphism of snow from 3-d images: Time evolution of microstructures, physical properties and their associated anisotropy. Cryosphere **8**, 2255–2274 (2014).

69. Wautier, A., Geindreau, C. & Flin, F. Linking snow microstructure to its macroscopic elastic stiffness tensor: a numerical homogeniza-tion method and its application to 3-D images from X-ray tomo-graphy. Geophys. Res. Lett. **42**, 8031–8041 (2015).

70. Haffar, I. et al. X-ray tomography for 3D analysis of ice particles in jet A-1 fuel. Powder Technol **384**, 200–210 (2021).

71. Fang, Q. & Boas, D. A. Tetrahedral mesh generation from volumetric binary and gray-scale images, 1142–1145 (IEEE Press, Boston, Mas-sachusetts, USA, 2009).

72. Cignoni, P. et al. Scarano, V., De Chiara, R. & Erra, U. (eds) Meshlab: an open-source mesh processing tool. (eds Scarano, V., De Chiara, R. & Erra, U.) Eurographics Italian Chapter Conference, 129–136 (The Eurographics Association, Salerno, Italy, 2008).

73. Dawson-Haggerty, M. trimesh, 3.2.0. https://trimsh.org/ (2019).

74. Kaempfer, T. U., Hopkins, M. A. & Perovich, D. K. A three- dimensional microstructure-based photon-tracking model of radiative transfer in snow. J. Geophys. Res. **112**, 24113 (2007).

75. Picard, G., Libois, Q. & Arnaud, L. Refinement of the ice absorption spectrum in the visible using radiance profile measurements in antarctic snow. Cryosphere **10**, 2655–2672 (2016).

76. Toon, O. B., McKay, C. P., Ackerman, T. P. & Santhanam, K. Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres. J. Geophys. Res. **94**, 16287–16301 (1989).

77. Carmagnola, C. M. et al. Snow spectral albedo at summit, Green- land: measurements and numerical simulations based on physical and chemical properties of the snowpack. Cryosphere **7**, 1139–1160, https://doi.org/10.5194/tc-7-1139-2013 (2013).

78. Martin, O. A., Kumar, R. & Lao, J. Bayesian Modeling and Computa-tion in Python (Boca Raton, Florida, USA, 2021).

79. Salvatier, J., Wiecki, T. V. & Fonnesbeck, C. Probabilistic programming in python using pymc3. PeerJ Comput. Sci. **2**, e55 (2016).

80. Hansen, J. E. & Travis, L. D. Light scattering in planetary atmo-spheres. Space Sci. Rev. **16**, 527–610 (1974).

81. Mishchenko, M. I. & Macke, A. Asymmetry parameters of the phase function for isolated and densely packed spherical particles with multiple internal inclusions in the geometric optics limit. J. Quant. Spectrosc. Radiat. Transfer **57**, 767–794 (1997).

82. Stamnes, K., Tsay, S. C., Jayaweera, K. & Wiscombe, W. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl. Opt. **27**, 2502–2509 (1988).

83. Gardner, A. S. & Sharp, M. J. A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization. J. Geophys. Res. **115**, https://doi.org/10.1029/ 2009jf001444 (2010).

84. Dang, C., Brandt, R. E. & Warren, S. G. Parameterizations for nar-rowband and broadband albedo of pure snow and snow containing mineral dust and black carbon. J. Geophys. Res.: Atmos. **120**, 5446–5468 (2015).

85. Hagenmuller, P. et al. Motion of dust particles in dry snow under temperature gradient metamorphism. Cryosphere **13**, 2345–2359 (2019).

86. Xiong, C. et al. A new hybrid snow light scattering model based on geometric optics theory and vector radiative transfer theory. IEEE Trans. Geosci. Remote Sens. **53**, 4862–4875 (2015).

87. Hagenmuller, P., Chambon, G., Flin, F., Morin, S. & Naaim, M. Snow as a granular material: assessment of a new grain segmentation algorithm. Granular Matter **16**, 421–432 (2014).

88. Hagenmuller, P., Matzl, M., Chambon, G. & Schneebeli, M. Sensi-tivity of snow density and specific surface area measured by microtomography to different image processing algorithms. Cryo-sphere **10**, 1039–1054 (2016).

Acknowledgements

This research has been supported by the Agence Nationale de la Recherche (MiMESis-3D project—grant no. ANR-19-CE01-0009). IGE and CNRM/CEN are part of Labex OSUG@2020 (Investissements d'Avenir - grant no. ANR-10-LABX-0056). Marie Dumont has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (IVORI - grant no. 949516). We are also grateful to the 3SR tomographic plateform and the ESRF ID19 beamline where all the 3D images used in this study have been acquired. Most of the computations presented in this paper were per-formed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. The authors would like to acknowledge Dr. Kenneth G. Libbrecht, Florent Domine and Météo-France for allowing the use of their photographs in this study.

Author contributions

A.R., G.P., and M.D. designed the study. G.P. and A.R. developed the ray-tracing model and the estimation method. F.F. conducted the X-ray tomography measurements, prepared the 3D mesh images, and obtained funding. A.R., G.P., M.D., and L.A. performed the analysis of the simulations. Q.L., G.P., A.R., and M.D. explored the relationship B = n². A.R. drafted the original manuscript, and all authors made suggestions for the improvement of the manuscript.

Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-39671-3.

**Correspondence** and requests for materials should be addressed to Alvaro Robledano.

**Peer review information** *Nature Communications* thanks Cenlin He, Petri Räisänen, Stephen Warren, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023