A generalized photon-tracking approach to simulate spectral snow albedo and transmissivity using X-ray microtomography and geometric optics

A majority of snow radiative transfer models (RTM) treat snow as a collection of idealized grains rather than an organized ice-air matrix. Here we present a generalized multi-layer photon-tracking RTM that simulates light transmissivity and reflectivity through snow based on X-ray microtomography, treating snow as a coherent structure rather than a collection of grains. Notably, the model uses a blended approach to expand ray-tracing techniques applied to sub-1 cm snow samples to snowpacks of arbitrary depths. While this framework has many potential applications, this study’s effort is focused on sim5 ulating light transmissivity through thin snowpacks as this is relevant for surface energy balance applications and sub-nivean hazard detection. We demonstrate that this framework capably reproduces many known optical properties of a snow surface, including the dependence of spectral reflectance on snow grain size and incident zenith angle and the surface bidirectional reflectance distribution function (BRDF). To evaluate how the model simulates transmissivity, we compare it against spectroradiometer measurements collected at a field site in east-central Vermont. In this experiment, painted panels were inserted at 10 various depths beneath the snow to emulate thin snow. The model compares remarkably well against the spectroradiometer measurements. Sensitivity simulations using this model indicate that snow transmissivity is greatest in the visible wavelengths and is limited to the top 5 cm of the snowpack for fine-grained snow, but can penetrate as deep as 8 cm for coarser grain snow. An evaluation of snow optical properties generated from a variety of snow samples suggests that coarse grained low density snow is most transmissive. 15

This model can be divided into two distinct components. The first of which determines key snow medium optical properties by launching photons into 3D closed-surface renderings of snow samples derived from µCT scans with a voxel resolution of ≈ 20 µm. The second uses the optical properties derived from the first part to drive a 1D photon-tracking model whereby individual photon packets are prescribed a random initial position and incident direction on the snow. Each individual photon packet then has a unique path whereby all of the energy contained within a given packet travels in the same direction and the amount of energy within a given packet is depleted continuously according to absorption within the medium. Note that for 95 both model components, the ice refractive indices reported by Warren and Brandt (2008) are used to compute scattering and absorption.

Snow Optical Properties
The 1D medium model requires three key optical properties: the extinction coefficient (γ ext ), the mean path fraction traveled within ice (F ice ), and the scattering phase function (p(cosΘ)). In considering light as a ray traveling through the snow medium, 100 which is scattered each time it intersects an air/ice boundary and partially absorbed within the ice, the extinction coefficient is related to the distance traveled between scattering and absorption events. The phase function determines the change in direction of the ray during a scattering event, and the ice-path fraction, when combined with the ice absorption coefficient, determines the mean energy depleted from the ray for a given distance traveled between scattering events. For a given snow sample, γ ext is determined following the method described in Xiong et al. (2015) applied to the 3D rendering of snow as opposed to 105 an idealized bicontinuous medium. In this framework, photons are initialized at a random position within the snow sample, and launched in a random direction for a specified distance (L). If the photon is initialized within the air, the probability of extinction (P ext ) is 1 if a boundary is intersected over L, otherwise it is 0. In the case where the photon is initialized within the ice medium, P ext =1 if a boundary is intersected over L, otherwise it is given as: 110 where κ λ is the wavelength-dependent absorption coefficient of ice, which is related to the imaginary part of the ice refractive index (k): This slight modification is made to account for the added probability of extinction due to absorption of the photon within the ice particle. Note that this also introduces a minor wavelength dependence into the extinction coefficient. Using this method, 115 a probability of extinction can be determined for distance L. This method is repeated for several distances ranging from the voxel resolution (20 µm) to the width of the snow sample volume (e.g., 10 mm). The extinction coefficient is then determined using a curve fit to Beer-Lambert law: 4 https://doi.org/10.5194/tc-2021-310 Preprint. Discussion started: 21 October 2021 c Author(s) 2021. CC BY 4.0 License.
The mean fractional ice path (F ice ) is determined by tracking individual photons as they travel throughout the aggregate 120 snow sample. This framework closely mimics that of Kaempfer et al. (2007) in that photons travel through the snow medium and change direction according to Snell's law of refraction and a probabilistic representation of Fresnel's law of reflectance.
Here, a photon is initialized at a random starting point somewhere within the snow sample and launched in a random direction.
The photon is tracked until it exits the medium, and the F ice is simply the ratio of the distance traveled within ice over the total distance traveled. This is repeated for a large number of photons to determine an average F ice .

125
Fresnel's law dictates that the fractional reflection and transmission of light at a boundary is related to the incident angle (θ i ) and the refractive indices (n) of the two media separated by the boundary: and

130
where R h and R v are the horizontally and vertically polarized reflectances. Assuming that the radiation is unpolarized (e.g., natural light), the reflectance (R) is: Through energy conservation, the transmittance (T ) is simply: 135 Then, if the vector normal to the boundary plane (v n ) is oriented towards the medium with refractive index n 1 , the direction unit vectors for transmitted and reflected radiation are computed as: and 140 wherev r andv t are the reflection and transmission unit direction vectors, respectively.
The phase function is determined by separating out individual snow grains from the reconstructed 3D snow sample and probing them with photons to compute the scattering angle where the cosine of the scattering angle (Θ) is the dot product between the directional unit vector of radiation incident on the particle (Ω ) and the directional unit vector of the scattered radiation (Ω) in the cartesian coordinate space: The phase function is constructed by first initializing a photon outside of a given particle and firing it in a random direction towards the particle. The photon then interacts with the snow particle, guided by Eqs. 4 -9. For each collision (i), the amount of energy exiting the particle into the directional bin (cosΘ j ) is tracked, and the remaining photon energy is depleted as energy exits, or is absorbed within, the ice particle. This is repeated for an arbitrary number of collisions (n) until there is less than 150 0.1% of the initial energy left (Fig. 1). For the initial collision (i = 0), the energy exiting the particle is simply the reflected fraction of the incident ray, and for each subsequent collision, it is the transmitted multiplied by the remaining ray energy (e.g., Malinka, 2014): where s n is the distance traveled within the particle between the boundary intersections n−1 and n. For visible wavelengths, 155 the impact of absorption on the phase function is negligible, however for wavelengths exceeding 1000 nm, it becomes more important. The resulting distribution of energy is converted to a phase function defined relative to the ray initially incident upon the scattering particle following Grundy et al. (2000): where N is the total photon energy and N j is the total photon energy directed into bin j.

160
In practice, while a majority of photons require only a few collisions to reach the 0.1 % energy threshold, we cap the number of internal collisions to a maximum of 10 to limit both computation time and the accumulation of error caused by imperfect Figure 1. Schematic illustration of the photon-tracking method used to determine the scattering phase function for a grain with an arbitrary shape.
3D particle representations. To evaluate this method, we use it to estimate the phase function of a collection of idealized crystal habits, including spheres, hexagonal plates, and columns with a size parameter of 1000 (Fig. 2). We find that these phase functions are in good agreement with Mie theory for spheres, and with phase functions reported for other shapes presented in 165 previous studies focused on the scattering properties of spherical, idealized, and irregular snow crystal shapes (Iaquinta et al., 1995;Macke et al., 1996;Yang and Liou, 1996;Grundy et al., 2000;Malinka, 2014;Ishimoto et al., 2018). To determine the phase function for 1D model, this ray-tracing method is applied to a selection of rendered grains from the µCT sample, and p(cosΘ j ) is computed for each bin following eq. 12 where N j is determined from all selected grains. Note that in this method the selected grains are checked to ensure that they contain no facets touching the sample boundaries.

1D photon-tracking Model
Once the required optical properties of the snow sample are determined by launching photons through µCT sample volumes, a 1D photon-tracking model is used to simulate snow spectral albedo, transmissivity, and Bidirectional Reflectance Distribution Function (BRDF). The 1D model is used in place of the explicit photon-tracking model described by Kaempfer et al. (2007) in order to allow for the computationally feasible simulation of spectral albedo and transmissivity for snow covers with 175 depths exceeding 1 cm with sufficient grain resolution. Additionally, it is used complications associated with lateral boundary treatment and stitching multiple µCT scans together into a coherent snow lattice. Our 1D model is based largely on Jacques (2010), which describes a semi-random 1D multi-layer Monte Carlo photon-tracking approach for application in the field of biomedical imaging. In this framework, discrete, plane-parallel, snow layers with optical properties constant throughout each layer are first prescribed. Then a photon packet is initialized at some starting position (X 0 ) with cartesian components of (x 0 ,  Phase functions for idealized snow particles following equations 11 and 12 for incident radiation at λ=900nm. The size parameter (kL) for each particle is set to 1000. Note that for the rendered grain, the approximate grain diameter yielded a size parameter of approximately 1400. 500,000 photons were used to generate the phase function. AR is the axis ratio between the long "c" and short "a" axes of the hexagonal column.
y 0 , z 0 ) and an initial energy of unity (E = 1). An initial unit direction vector for the photon is given in cartesian coordinates as: where θ is the solar zenith angle, and φ is the azimuth angle clockwise from x. This initial direction can be prescribed randomly (i.e., diffuse radiation), or at any specified zenith/azimuth angle (i.e., direct radiation), or as a mixture of both diffuse and direct 185 radiation.
Once the initial position is set, the photon is launched into the medium, and travels a distance s before experiencing a scattering event. s is computed statistically using the Beer-Lambert law and the medium extinction coefficient (Jacques, 2010): where ζ is a random uniform number between 0 and 1. The new position in the medium is: At the scattering event, the photon packet is given a new direction unit vector according to the scattering phase function.
Because this framework treats the scattering phase function as a probability distribution function (PDF), the scattering angle Θ is determined by choosing a random sample from p(cosΘ) PDF:

195
where P is the probability of light being scattered into a cone with solid angle dΩ in the direction Θ from the incident radiation given the phase function.
Then the new direction vector is determined from Θ (Jacques, 2010): where φ is given as a uniform random number between 0 and 2π, the 0 subscript represents the incident direction, and µ x , µ y , 200 and µ z make up the components of the unit direction vector.
Photon energy is depleted over distance s according to the ice absorption coefficient and F ice as determined from the µCT data instead of using a medium absorption coefficient: where E is the new photon energy, and E 0 is the incident photon energy.

205
To achieve statistical energy conservation, a "Russian Roulette" function is used to determine whether or not to fully absorb (i.e., kill) the photon packet once its energy falls below a prescribed threshold (Iwabuchi, 2006;Jacques, 2010). This is given as: where ζ is a random number between 0 and 1, and m is a prescribed constant on the order of 1-10. By treating absorption continuously rather than probabilistically, the number of photons required to attain a robust solution is significantly reduced, and further ensures that the model cannot get stuck in an infinite loop.
If the z position of a photon-packet is above the top of the snow surface (i.e., it has exited the top of the snowpack), the remaining energy within the packet is added to the total reflected energy and the photon is eliminated. In an open lowerboundary configuration, if a photon-packet z position is less than 0 (i.e., it has exited the bottom of the snowpack) the remaining 215 energy is added to the total transmitted energy, and the photon is eliminated. Alternatively, a lower boundary can be simulated with a specified spectral reflectance such that a portion of the photon energy will be absorbed at the lower boundary, and the remaining energy will be reflected upward. Once all photons have been eliminated from the model, the simulation is complete.
This model is extended to a multilayer configuration, by simply defining unique optical properties corresponding to specified depths throughout the snowpack. When a photon packet travels from one layer to another, its trajectory and energy depletion 220 are determined by the optical properties of the new layer.
The basic premise of this model is illustrated in figure 3, which traces the position and energy of two photons on a 2D plane as they travel throughout an idealized two-layer snowpack 10 cm deep.

Directional Conic Reflectance Function
The reflectance of a surface is often described using the concept of a BRDF (e.g., Stamnes and Stamnes, 2016). This concept 225 essentially represents a PDF of reflected direction of a ray of light impacting the surface from a given incident direction, and is used to simplify the complex reflectance properties of a rough surface (e.g., shading and multiple reflections). To estimate the BRDF from this model, we follow the methods described in Kaempfer et al. (2007). In this framework, the BRDF for specified incident zenith and azimuth directions is approximated using the Directional Conic Reflectance Function (DCRF), which computes the energy reflected into a cone in the direction: θ r ,φ r subtended by solid angle: dΩ: where I is the radiative flux, and the subscripts i, and r correspond to the incident and reflected radiation, respectively.

Snow sampling and spectroradiometer measurements in the field
To evaluate the model, we collected snow samples and snow surface spectroradiometer measurements at Union Village Dam (UVD) in Thetford, Vermont several times throughout the 2020-21 winter. The UVD site is a broad flat clearing surrounded by 235 deciduous forests spanning approximately 40000 m 2 , and bounded on the southern end by the Ompompanoosuc River. During each data collection, a snow pit was excavated and standard snow characteristics, such as snow depth, density, and grain size were measured manually. Several snow samples were carefully extracted in columns adjacent to the snow pit sidewalls in cylindrical containers 7 cm high x 1.9 cm in diameter (Fig. 4). These samples were transported in a hard, plastic cooler for 10 miles from the UVD site to the Cold Regions Research and Engineering Laboratory (CRREL). At CRREL they were stored in 240 a -20 • C cold room prior to µCT analysis. These samples were not casted (i.e. not preserved using a pore-filler).
Spectral reflectance and transmissivity data were collected using a Malvern Panalytical ASD FieldSpec 4 Hi-Res: High Resolution Spectroradiometer. The FieldSpec 4 has a spectral range of 350-2500 nm and a spectral resolution of 3 nm in the visible and 10 nm in the SWIR. The data collection was performed within 1.5 hours of solar noon in order to limit high zenith angle impacts. An optimization was conducted prior to the start of data collection and any time lighting conditions changed 245 in order to ensure accurate reflectance readings. Data collections were taken 2.5 to 3 feet above the snow surface using a 5 degree field of view optic lens, resulting in a measurement footprint diameter of approximately 6 cm. The collection strategy employed included taking a white reference reading from a pure reflective panel and five readings at different locations on the target surface; the mean of the five readings was used as the reflectance value for that specific location.
In this paper, we focus specifically on data collected on 12 February, 2021 as this day had the most stable ambient lighting 250 conditions and resulted in the majority of our snow and reflectance measurements. At the time of the measurements the sky was covered with a high optically thick overcast, and as a result the ambient lighting conditions were generally diffuse. The gray-scale horizontal slices using NRecon software (Bruker), which utilizes a modified Feldkamp cone-beam algorithm to produce a vertical stack of gray-scale cross-section images. Image reconstruction processing included sample-specific post alignment, Gaussian smoothing using a kernel size of 2 to reduce noise, sample-specific ring artifact correction of dead pixels, 270 beam hardening correction, and X-ray source thermal drift correction. A cylindrical volume of interest with a diameter of 1.6 cm was selected from the scanned samples in order to eliminate edge effects caused by the sampling process.
Resulting grayscale images are segmented into two phases: air (lowest X-ray absorption), and snow (highest X-ray absorption). Segmenting thresholds for each phase are determined by finding the local minimum between peaks on the histogram showing all grayscale values, and using that value as a global threshold for each scanned sample. The resulting binarized data 275 are despeckled so that any objects less than 2 pixels in diameter were removed.
The final binarized images are then used to construct 3D representations of dry snow samples for input into the RTM. This  To build a full sample mesh, the grains are processed on an individual basis and then combined together to create the full mesh. To process each individual grain, a contour-based surface reconstruction process was developed to generate grain 290 surfaces from the voxels that make up the grain. This method uses a subset of the binary sample array that contains the target grain, including both snow and adjacent air voxels. The subset array is then refined to increase the resolution. A Gaussian filter is applied to smooth the refined array, diminishing pixelated appearance of the voxelized snow-air interface, producing a smooth level set from which to extract the grain surface (Fig. 6). The smoothed level set is then used to define an isosurface at the snow-air boundary, providing control over where the boundary is drawn with respect to the voxels.

295
Finally, to extract the isosurface from the 3D voxel array, we apply the Marching Cubes method. In this technique, the input volume is divided into a discrete set of cubes. The algorithm then determines how the surface intersects with a given cube,  1994; Natarajan, 1994;Chernyaev, 1995;Lewiner et al., 2003). For this work, we used the adaptation implemented by Lewiner et al. (2003), which improved the algorithm to resolve face and internal ambiguities, extended the lookup table, and guaranteed correct topology. As a final step, each grain is "repaired" to remove any defects and degenerate elements and ensure a manifold 305 surface according to Attene (2010), and then decimated to reduce the overall number of triangles that comprise the surface thereby lowering the computational requirements. Overall, this method appears to accurately characterize the snow within the µCT sample with computed mesh snow sample densities within 1.5% of snow densities computed from the raw voxels. Figure   7 shows a 2D cross section comparing grain boundaries to the raw pixels of the image and selected example 3D rendered grains are show in Figure 8.

General Evaluation
An initial evaluation of the model is performed by simulating the spectral albedo for two idealized 20 cm deep snowpacks with uniform optical properties throughout. For these snowpacks, the optical properties are determined from 3D meshes generated by two characteristically distinct µCT samples. One mesh is representative of fresh, fine-grained snow near the surface, and 315 the other of large facets near the bottom of the snowpack (Fig. 9). For each mesh, the total mesh volume is approximately    Table 1. For each sample, the spectral albedo is computed for wavelengths between 400-1300 nm at 25 nm intervals with diffuse incident radiation. This comparison demonstrates that the model capably reproduces some known behavior of spectral albedo, and specifically, it shows a strong wavelength dependence on snow microstructure that favors the NIR (Fig. 10a). The spectral albedo is relatively uniform 320 between the two snowpacks for the spectral range between 400 and 800 nm, and then the albedos diverge, with a more rapid decrease in albedo for the coarser-grained snow.
We then assess the relationship between simulated spectral albedo and incident zenith angle for the fine grain snow sample at four different wavelengths to evaluate the model's ability to simulate anisotropy in the surface reflectance (Fig. 10b). This analysis shows an exponential increase in albedo at high zenith angles that is most pronounced in the NIR, consistent with 325 observed behavior. This result indicates that the model is capable of simulating surface anisotropy with good fidelity. As a related evaluation, the model-simulated DCRF is computed as a function of zenith angle (Fig.11). This analysis reveals that the reflectance is mostly isotropic for zenith angles less than approximately 55 • at which point the surface becomes increasingly forward scattering, consistent with previous observational and modeling studies (Kaempfer et al., 2007;Dumont et al., 2010;Xiong et al., 2015;Jiao et al., 2019).

330
Finally, we use the model to provide an initial assessment of the impacts of snow microstructure on simulated spectral transmissivity. To accomplish this, the optical properties of the µCT samples in Fig. 9 are used to simulate and compare the spectral transmissivity at varying depths (Fig. 12). The transmissivity is highest at the short, non-absorptive, wavelengths and gradually decreases throughout the NIR, broadly matching quantitative snow transmissivity results reported in Perovich (2007) and Libois et al. (2013). The depth of the 5% transmissivity contour for the fine-grain snow sample is approximately 2.5 cm for 335 the visible, and decreases to approximately 1.5 cm for the NIR (Fig. 12a), indicating that the fine-grain snow optical thickness is on the order of only a few centimeters. In contrast, the transmissivity for the coarse grain snow is increased near the surface, and the depth of the 5 % contour increases to 7.5 cm for the visible and 3 cm for the NIR (Fig. 12b).

Evaluation against UVD Data
The optical properties used in the 1D RTM were determined from four approximately 800 mm 3 µCT samples, with each 340 sample representing a 2 cm thick layer within the top 8 cm of the snowpack. The RTM is then configured with 4 layers according to these optical properties (given in Table 2). The top three layers are each 2 cm thick, and the bottom layer is 28 cm thick, such that the entire snow depth amounted to 34 cm. We choose this configuration working under the hypothesis that the snow microstructure below 8 cm had little impact on the measured surface spectral albedo. To simulate the panels, the snowpack depth is modified to be 4.75 and 2.5 cm deep with a 100 % absorptive lower boundary while maintaining the 345 layering corresponding to Table 2.
There is remarkably good agreement between our observations and the model (Fig. 13) and in particular, the model accurately simulates the impact of the inserted panel on the surface albedo for wavelengths shorter than 1000 nm for both the 4.75 and 2.5 cm depths. The model and observations diverge after 1200 nm, in particular, the simulated albedo is substantially higher in this range than measured. We hypothesize that this is primarily due to the phase function approximation computed from the rendered 350 snow grains. Additional simulations that use the phase function computed from idealized spheres and columns support this difference between the coarse grain and fine grain samples as a function of wavelength and depth. The 5% transmissivity contour is plotted for the fine grain (black) and coarse grain (red) samples.  γ ext and F ice . This analysis is performed using several µCT sample volumes collected on different dates, at different locations, 360 and for various snow types. Note that each µCT sample is approximately 800 mm 3 and the sample SSA and ρ s are determined from the µCT 3D rendering. This analysis reveals that F ice has a very robust relationship with snow density (Fig. 14a) described by the linear fit: with r 2 = 0.81. In contrast, while γ ext generally seems to increase as a function of SSA, this relationship is not as well 365 constrained as the relationship between ρ s and F ice (Fig. 14b). Instead, we find that γ ext is generally better approximated by a multivariate regression function that includes both SSA and density: and increases r 2 to from 0.25 to 0.79. This relationship dictates that the extinction coefficient increases both as a function of SSA and ρ s , indicating that the highest extinction coefficients will be associated with small, tightly packed snow grains. This is qualitatively consistent, though not directly comparable, with analytical formulations in Kokhanovsky and Zege (2004), who show the extinction coefficient is related to both surface area and grain concentration, which is a good proxy for snow density.

375
In contrast, with ice-crusts the extinction coefficient is more strongly related to SSA, and largely independent of ρ s , leading to low extinction coefficients despite very high snow densities, counteracting the relationships described in equations 22,23. This behavior can be attributed to the fact that for crust layers, snow density is no longer a good proxy for particle concentration and therefore does not significantly affect the extinction coefficient. Overall, these results suggest that commonly observed physical snow properties can be used to approximate optical properties in conjunction with an appropriate phase function for 380 the medium photon-tracking RTM model for non-crust snow layers.
To further assess how these two specific snow optical properties, γ ext and F ice , affect the greater simulated spectral transmissivity, we perform a sensitivity analysis by comparing the 5% transmissivity contour depth for three fractional ice paths: 0.30, 0.46, 0.69 at two fixed γ ext values: 2.34, 0.81 mm −1 . The two γ ext values correspond to the max, min values found in the previous analysis and presented in Fig. 14b. The three F ice values correspond to the max, min, and mean values (Fig. 14a). We 385 compare the influence of F ice at both the max and min γ ext values, since we anticipate the strength of its influence will vary according to γ ext . While we pair the maximum γ ext with the F ice , we note that high values ofγ ext are more likely to coincide with high values of F ice due to the shared dependence of these variables on snow density in most snowpacks.
The results of this analysis indicate that γ ext has a much larger impact on snow transmissivity than F ice . This is unsurprising due to the highly scattering nature of snow. In particular the snow medium is much more transmissive at the minimum γ ext 390 value with the 5 % transmissivity contours exceeding 7 cm in the visible and 3 cm in the NIR. This is in contrast to the transmissivity of the maximum extinction coefficient, which never exceeds 4 cm of depth (Fig. 15). F ice can either amplify or dampen the effect of γ ext in the NIR, with lower values of F ice leading to in increase in transmissivity at a given depth compared to high values of F ice .

395
One key objective of this work is to expand beyond determining snow optical properties from a specified distribution of grains with idealized shapes, and instead represent snow as an organized structure in determining them. One snow type where this approach may be particularly advantageous is in understanding the optical properties of snow crust layers, which do not fit easily into the collection of particles or air bubbles approximation. This is supported by our finding that γ ext is very well highest snow densities and the lowest extinction coefficients and SSAs. Another snow type this approach may be well suited for is highly aged snow that has metamorphosed into a collection of large irregularly shaped grain clusters and pore spaces.
While ideally, the optical properties would vary slightly with wavelength due to the impact of absorption on γ ext and the phase function, we chose to leave the optical properties independent of wavelength. This choice is made entirely to reduce the computational burden of running the photon-tracking model for several wavelengths. Cursory sensitivity tests performed 405 to assess the impact of this choice on the optical properties supported the use of wavelength-independent optical properties, as both p(cosΘ) and γ ext exhibited generally a negligible dependence on wavelength for λ > 1400. We suspect that this is due to the fact that a ray is much more likely to be scattered at an air/ice boundary than by absorption within the particle (i.e., snow has a high single scattering albedo). This may not be the case for all snow types, in particular for the NIR wavelengths within very large grains, and is worthy of future exploration. In this work we have presented a blended photon-tracking radiative transfer model in an effort to better understand the complicated influence of snowpack microstructure on snow spectral transmissivity in the geometric optics limit. A foundational goal of this modeling approach is to expand upon previous approaches aimed at incorporating 3D renderings of real snow microstructure into radiative transfer models for snowpacks of arbitrary depth, while maintaining the Monte Carlo aspects of 415 the model. To accomplish this, existing methods for simulated photon interactions with rendered elements are employed to determine key optical properties of the snow (Grundy et al., 2000;Kaempfer et al., 2007;Xiong et al., 2015).
An evaluation of this framework for consistency with known behavior of spectral snow albedo revealed that this framework can successfully reproduce the dependency of spectral albedo and grain size, as well as the surface anisotropy at high incident zenith angles. Furthermore, a comparison of the simulated snow albedo against spectroradiometer measurements collected in 420 the field over snow with varying depths indicate that the model can simulate the effects of an underlying surface on spectral albedo with high accuracy.
In comparing two different snow samples, it was revealed that snow microstructure has a large impact on snow transmissivity in the visible spectrum and near the snow surface, increasing the 5 % transmissivity depth from approximately 4 cm for a fine grain snow sample to 7.5 cm depth for a coarse grain sample. A brief sensitivity analysis of the optical properties revealed that 425 lowering the medium extinction coefficient acted to reduce the albedo and increase transmissivity in the visible bands, while the fractional ice path (F ice ) impacted the rate at which albedo and transmissivity decreased as a function of wavelength in the NIR. Accordingly, we anticipate that snowpacks made up of large grains with low fractional ice-paths will be the most transmissive.
Overall, while current efforts are focused on using this model to better understand snow transmissivity, it shows promise as 430 a broadly applicable snow RTM that has a strong direct connection to µCT snow samples. While currently, it is limited to the geometric optics approximation for clean snow and unpolarized radiation, ongoing and anticipated future efforts are aimed at improving the grain segmentation and rendering process, incorporating polarization, parameterizing diffraction, and including light absorbing particulates (LAPs). In particular, recent multiphase image segmentation techniques (e.g., West et al., 2018) could be used to better separate snow, air, and LAPs in a µCT sample allowing for the impact of LAPs to be determined through 435 ray-tracing. Furthermore, because the model operates entirely as a photon-tracking model, it is a natural fit with macro-scale ray-tracing and therefore could be used to investigate the reflectance of rough snow surfaces such as sun cups or sastrugi.
Code and data availability. The mesh generation and RTM code with associated documentation is available in preliminary 'as is' format on Github at (lhttps://github.com/wxted/CRREL-GOSRT.git). Sample data files used to generate figure 9a is available on Github as sample data. Additional limited sample data, including rendered microCT meshes, and spectroradiometer data used for this paper are available upon Tech. Note TH-387+ STR, 1993.