Mapping surface hoar from near-infrared texture in a laboratory

Surface hoar crystals are snow grains that form when water vapor deposits on the snow surface. Once buried, surface hoar creates a weak layer in the snowpack that can later cause large avalanches to occur. The formation and persistence of surface hoar are highly spatiotemporally variable, making its detection difficult. Remote-sensing technology capable of detecting the presence and spatial distribution of surface hoar would be beneficial for avalanche forecasting, but this capability has yet to be developed. Here, we hypothesize that near-infrared (NIR) texture, defined as the spatial variability of reflectance magnitude, may produce an optical signature unique to surface hoar due to the distinct shape and orientation of the grains. We tested this hypothesis by performing reflectance experiments in a controlled cold laboratory environment to evaluate the potential and accuracy of surface hoar mapping from NIR texture using a near-infrared hyperspectral imager (NIR-HSI) and a lidar operating at 1064 nm. We analyzed 41 snow samples, three of which were surface hoar and 38 of which consisted of other grain morphologies. When using NIR-HSI under direct and diffuse illumination, we found that surface hoar displayed higher NIR texture relative to all other grain shapes across numerous spectral bands and a wide range of spatial resolutions (0.5–50 mm). Due to the large number of spectral- and spatial-resolution combinations, we conducted a detailed samplewise case study at 1324 nm spectral and 10 mm spatial resolution. The case study resulted in the median texture of surface hoar being 1.3 to 8.6 times greater than that of the 38 other samples under direct and diffuse illumination (p < 0.05 in all cases). Using lidar, surface hoar also exhibited significantly increased NIR texture in 30 out of 38 samples, but only at select (5–25 mm) spatial resolutions. Leveraging these results, we propose a simple binary classification algorithm to map the extent of surface hoar on a pixelwise basis using both the NIR-HSI and lidar instruments. The NIR-HSI under direct and diffuse illumination performed best, with a median accuracy of 96.91 % and 97.37 %, respectively. Conversely, the median classification accuracy achieved with lidar was only 66.99 %. Further, to assess the repeatability of our method and demonstrate its mapping capacity, we ran the algorithm on a new sample with mixed microstructures, with an accuracy of 99.61 % and 96.15 % achieved using NIR-HSI under direct and diffuse illumination, respectively. As NIR-HSI detectors become increasingly available, our findings demonstrate the potential of a new tool for avalanche forecasters to remotely assess the spatiotemporal variability of surface hoar, which would improve avalanche forecasts and potentially save lives.

1 Introduction and background

Mountainous and polar snowpacks are commonly blanketed by surface hoar, unique ice crystals that grow when water vapor deposits on the snow surface (Horton and Jamieson, 2017). Hoar crystals can grow to several centimeters in length and stand in a predominately vertical (but often quite variable) orientation atop the snow surface. Surface hoar is a prominent concern for avalanche forecasters due to its propensity for creating weak layers in the snowpack that can cause large avalanches. Generally, a weak layer is formed due to a snow layer being weakly bonded to the slab above, as in the case of an ice lens, or being of low shear strength, such as with surface hoar. Once buried, surface hoar layers are prone to fracture propagation and avalanche release (Horton and Jamieson, 2017; Jamieson and Schweizer, 2000). For instance, Birkeland (1998) found that nearly one-third of the large natural avalanches in southwestern Montana failed on buried layers of surface hoar, while Jamieson and Schweizer (2000) similarly observed in a Swiss dataset that 40 % of avalanches released on a weak surface hoar layer.

Surface hoar formation and persistence is highly spatiotemporally variable. Its presence is difficult to predict because it depends on the complex spatial distributions of the atmospheric water vapor, precipitation, wind, radiation, and vegetation that are present in polar and mountainous environments. Temporally, surface hoar can be promptly destroyed by environmental influences, such as wind, or it can persist for weeks (Champollion et al., 2013; Lutz and Birkeland, 2011). The formation of surface hoar can occur across the landscape or in isolated sub-slope regions. Therefore, the detection of surface hoar before it becomes a buried weak layer in the snowpack is an ideal application for remote sensing.

Optical remote-sensing retrievals of snow and ice commonly use near-infrared (NIR) wavelengths, where ice is absorptive and reflectance is sensitive to microstructure. In the NIR, snow grain size is the primary driver of snow reflectance and albedo; the increased path length of light through ice yields a reduction in reflectance. Despite the complex shapes of snow grains, snow is commonly represented as a collection of ice spheres with radius r_e, known as the optical or effective snow grain size, when simulating snow reflectance and albedo (Grenfell and Warren, 1999). The effective grain size is related to the physical snow microstructure through the ice surface-area-to-volume ratio (or specific surface area; SSA) by the following relationship:

$$\begin{array}{}\text{(1)}& r_{\mathrm{e}}={\displaystyle \frac{\mathrm{3}}{\mathrm{SSA}}}.\end{array}$$

In practice, the non-linear inverse relationship between NIR reflectance and r_e is leveraged to map r_e from reflectance measurements (e.g., Nolin and Dozier, 2000). Because grain size controls the broadband NIR albedo, estimates of r_e have historically dominated physical snow surface characterization, constituting a major goal of snow optics.

In recent decades, the use of NIR hyperspectral imaging (NIR-HSI) and (to a lesser extent) light detection and ranging (lidar) has enhanced r_e-mapping efforts. Hyperspectral instruments have superior spectral resolution relative to broadband or multispectral reflectance measurements, producing a nearly continuous measured spectrum for each pixel in an image. Lidar, on the other hand, holds key advantages as an active remote sensor. Lidar units scan a surface by emitting rapid pulses of light (most commonly at an NIR wavelength) and record both the relative strength of backscattered light that has reflected off a target and the two-way travel time. These two measurement types can be used independently; the travel time is used to measure the surface topography, whereas the backscattered magnitude has been demonstrated to be a useful measure of optical properties. Although far less validated than traditional passive reflectance measurements, lidar backscatter may provide adequate r_e estimates (Yang et al., 2017).

Despite advances in snow optics, Horton and Jamieson (2017) note the fundamental disconnect between snow surface characterizations conducted by the remote-sensing community (i.e., r_e mapping), and the physical properties relevant to avalanche release. For avalanche forecasters, characterizing snow surface microstructure with the morphological grain shapes defined in the International Classification for Seasonal Snow on the Ground (ICSSG; Fierz et al., 2009) and their associated mechanical properties is critically important because these mechanical properties help determine how well new snow will bond to the old snow surface. Hence, a forecaster would much prefer a map identifying a potential future weak layer, like surface hoar, than a map of r_e. Relating morphological grain shape to r_e is difficult because the effective grain size only considers the path length of ice, which is a complex function of grain shape, traditional grain size, bulk density, and other physical characteristics. Therefore, r_e is not particularly useful for avalanche forecasting operations, although room for this adaptation does exist. While certain magnitudes of r_e are generally related to grain shape (Domine et al., 2007; Matzl and Schneebeli, 2010), few studies have formally attempted to use NIR reflectance for mapping the morphological grain shape instead of r_e. Further, the studies that have attempted this (e.g., Bühler et al., 2014; Horton and Jamieson, 2017) have found that surface hoar crystals produce moderate reflectance signatures relative to other grain shapes, making them difficult to delineate from less-concerning snow microstructures based on NIR reflectance magnitude.

To the best of our knowledge, the only study to date that successfully identified surface hoar formation from NIR remote sensing was conducted by Champollion et al. (2013) in Antarctica. As opposed to evaluating magnitudes of NIR reflectance, the researchers leveraged an NIR texture signature, defined as the localized spatial variability in reflectance, to classify the presence of surface hoar using an infrared camera and an 850 nm artificial light source. Simply put, the researchers found that a large, localized variance in NIR reflectance, as quantified by a contrast index, was strongly correlated with surface hoar crystals. Similarly, in the preliminary findings presented by Walter et al. (2023), the researchers quantified a 600 % increase in NIR reflectance spatial variability during surface hoar formation, as measured with a 905 nm lidar unit.

Figure 1 Profile view of surface hoar atop an underlying snow layer (a) juxtaposed with an idealized schematic of light scattering when encountering a surface hoar layer (b). Because surface hoar crystals are typically large and vertically oriented, incoming photons may experience a very long path length and thus substantial absorption before reaching the underlying snow layer (leftmost case). Depending on the angle of interaction, that path may be considerably shortened (middle left case), or, as these crystals tend to be modestly spaced, photons may evade the surface hoar crystals entirely and pass straight into the underlying snow layer (middle case). Further, specular contributions are thought to increase in surface hoar layers (rightmost cases), which can produce particularly high or low reflectance depending on the angle of interaction. Last, because the surface hoar crystals will almost certainly have a different SSA than the underlying layer, light scattering back to the sensor directly from the surface hoar layer should further enhance the spatial variability of reflectance.

Here, we postulate that the physical phenomena for these findings include increasing specular contributions with surface hoar growth (Walter et al., 2023) as well as variable ice absorption and a variable path length (Fig. 1b). Depending on how an incoming photon interacts with the relatively large, often vertically oriented surface hoar crystals, the photon could experience a wide range of ice path lengths and thus absorption. Further, it is thought that surface hoar plates promote specular reflectance, which can cause large portions of radiation to either return to the sensor or to scatter forward toward the snowpack and into the “radiation trap” of hoar crystals. As a result, we hypothesize that the presence of surface hoar will coincide with a quantifiable increase in localized reflectance variance, a texture signature which could be used to map its distribution. If this is the case, it would enable fine spatial- and temporal-resolution mapping of surface hoar extent (prior to burial) in challenging-to-observe environments, particularly as NIR remote sensors and uncrewed aerial vehicles (UAVs) become more cost-effective as forecasting tools. However, such a texture analysis has never been fully evaluated, rigorously quantified for accuracy, or compared to a wide variety of well-defined microstructures to ensure that the NIR texture is indeed a unique defining feature of surface hoar.

To address this knowledge gap, we performed controlled cold laboratory experiments to determine whether NIR texture can delineate the extent of surface hoar. We first created snow samples with varying grain shapes and physical properties and quantified their microstructures using X-ray computed microtomography (micro-CT). Using both a compact NIR-HSI and a terrestrial lidar unit, we scanned each sample under a variety of illumination conditions. We subsequently analyzed the resulting maps of reflectance to produce measurements of NIR texture. Finally, we assessed the statistical significance of increased texture in surface hoar samples. We used our results to inform optimal thresholds for classifying surface hoar on a per-pixel basis before analyzing the accuracy of our resulting classified data products.

2 Methodology

We aimed to prepare laboratory snow samples with a wide variety of well-defined grain shapes and microstructures, acquire optical measurements, and perform a texture analysis to allow the delineation of surface hoar from other snow surface grain shapes. Section 2.1 describes snow sample preparation and physical characterization, Sect. 2.2 describes the acquisition of NIR-HSI and lidar data, Sect. 2.3 covers image texture analysis, Sect. 2.4 describes classification, and Sect. 2.5 evaluates the repeatability of our work.

Figure 2 Microscopy images of grains from each initial batch (left columns) and binary micro-CT cross-sections from representative samples (right columns). In the microscopy images, the grid size on the underlying blue grain card is 2 mm.

Table 1 Physical snow sample characteristics organized by primary grain shape and listed in order of decreasing surface-area-to-volume ratio. Grain shape abbreviations follow the ICSSG (Fierz et al., 2009).

2.1 Sample preparation and physical characterization

2.1.1 Sample preparation

We utilized Montana State University's Subzero Research Laboratory (SRL), a controlled cold laboratory environment, for sample preparation, testing, and assessment. The snow used in these experiments was a combination of laboratory-grown crystals produced in the SRL's snowmaking apparatus, which is similar to the systems presented in Schleef et al. (2014) and Abe and Kosugi (2019), and natural undisturbed snow that we collected from the surrounding area. To ensure the snow was completely dry, we kept all samples in a cold room at −30 °C and allowed them to equilibrate for at least 24 h prior to evaluation. We prepared 41 snow samples from 12 batches of differing snow grains (Fig. 2). From the bulk batches, we sieved snow grains through differing mesh sizes to further promote disparate microstructures (Table 1). The exception to this was surface hoar, which was grown in the laboratory atop rounded grains following the methods used by Stanton et al. (2016) and employing an apparatus similar to that used by Chandel et al. (2023) and Ozeki et al. (2020). Samples 25 and 26 consisted of in situ surface hoar growth at differing stages and measuring ∼ 0.5–2.0 cm in length. Meanwhile, Sample 24 featured these surface hoar grains redistributed through a large (6.30 mm) sieve in an effort to further examine the optical behavior of these grains after disrupting their typical vertical structure. Sample grain shapes included precipitation particles (PP), decomposing and fragmented precipitation particles (DF), rounded grains (RG), melt forms (MF), faceted crystals (FC), depth hoar (DH), and surface hoar (SH) (Fierz et al., 2009). We prepared snow samples that were microstructurally homogeneous, both laterally across the sample and vertically over the sample depth, in a rectangular 38 cm × 23 cm sample holder. The sample thickness was 3.8 cm, over twice the largest optically active depth estimates for our sample microstructures in the NIR spectrum (Nolin and Dozier, 2000).

2.1.2 Physical characterization

We thoroughly characterized the physical properties of each sample, as summarized in Table 1. First, we performed microscopy on representative grains from each batch prior to sieving (Fig. 2) and classified grain shapes using a crystal card and lens following Fierz et al. (2009). After sieving and sample preparation, we collected micro-CT data from each sample using a Bruker SkyScan 1173 housed in a −10 °C chamber within the SRL, generally following the protocol outlined by Donahue et al. (2021). To prepare samples, we used a cylindrical holder with 3 cm diameter × 4 cm length, which allowed for a voxel size of 14.5 µm. The voxel size of 14.5 µm was the finest spatial resolution achievable with the relatively large cylindrical micro-CT sample holder used in this study. The larger micro-CT sample holder was chosen to provide sufficient surface area for larger-grained samples, namely surface hoar, to be encapsulated and transported to micro-CT for measurement. We intuitively note that this will have resulted in the neglect or mischaracterization of fine dendrites smaller than this size. Consequentially, it is possible that the SSA of samples with a PP primary grain habit (particularly Samples 1–5) was underestimated by micro-CT, although this is of little consequence for the optical texture analysis presented here. We obtained measurements using a 42 kV, 190 mA X-ray beam with 100 ms exposure time, with each sample rotated 180° in 0.7° increments. After scanning, we performed thresholding of greyscale images into ice and air phases by visual inspection (e.g., Fig. 2) and used a despeckling filter to remove white and black speckles. Reconstructions performed via the marching cubes method (Lorensen and Cline, 1987) allowed us to determine the volume and surface area in 3D, which we used to calculate the SSA and density of each sample.

Figure 3 Laboratory data collection schematic for hyperspectral imaging under direct (a) and diffuse (b) illumination conditions as well as for lidar (c). Data regions of interest (ROIs) within the snow sample are illustrated in panel (d). Angle Θ refers to the viewing incidence angle (as well as the illumination incidence angle when both the lidar and NIR-HSI are under direct illumination).

2.2 Optical measurements

We examined NIR texture under a variety of illumination and viewing conditions using both lidar and NIR-HSI independently. For NIR-HSI, we constructed laboratory setups for both direct (see the “Direct-illumination data collection” section) and diffuse (see the “Diffuse-illumination data collection” section) illumination conditions. As lidar (Sect. 2.2.2) produces its own direct illumination via laser irradiance, this was the only possible illumination condition for lidar analysis. Furthermore, it is well-understood that the incidence angle impacts snow reflectance in the NIR spectral region; snow is nearly Lambertian under nadir incidence and predominantly forward scattering at off-nadir incidence. To consider this, for each of the aforementioned illumination configurations/instruments, we collected data with the sample container (1) perpendicular to the detector (nadir) and (2) tilted 10° away from the detector. Hereafter, these incidence angles of 0 and 10° are termed Θ. The result was six datasets representing each illumination–instrument–incidence combination for every sample. We acquired all optical data immediately prior to micro-CT analysis at a constant temperature of −10 °C.

2.2.1 Hyperspectral imaging

We used a Resonon Inc. Pika NIR-640 near-infrared hyperspectral imager to map snow reflectance in the NIR (http://www.resonon.com, last access: 11 November 2023). Donahue et al. (2021) provide a detailed description of the instrument. Briefly, the spectral resolution of the imager ranges from 2.39 to 2.50 nm and it measures 336 bands across the NIR region from 891–1711 nm. It constructs a 2D image containing the full spectrum in each pixel by collecting the image line by line, known commonly as a “push broom” or “line” scanner. Thus, to collect an image, the camera must move (translating or rotating) relative to the scene or the scene must move relative to the imager. We used a Resonon benchtop linear scanning stage to move the sample beneath the sensor.

Direct-illumination data collection

We positioned the hyperspectral imager above the linear translating stage that held the samples. For more details on the benchtop apparatus, see Donahue et al. (2022). The lens of the imager is surrounded by a set of four halogen lamps that produce direct illumination (Fig. 3a). The halogen lamps and lens of the imager are at a height of 38 and 47 cm above the snow surface, respectively. We used a large Spectralon white diffuse-reflectance panel to perform calibration, resulting in the measured reflectance factor (R) for each band in every individual pixel of the image. The Spectralon panel is 30.5 cm × 30.5 cm and thus larger in both dimensions than our optical ROI. We built a sample holder that had the same external dimensions as our snow sample holders but was specifically made to hold the Spectralon panel such that it was both centered on the ROI and at the same distance from the illumination source as the snow surfaces. For each snow sample scan, we also conducted a reference scan with the Spectralon panel. This allowed for pixel-by-pixel calibration of the entire optical ROI, thus accounting for any heterogeneous illumination. We made these reference measurements for each sample and each illumination condition.

Diffuse-illumination data collection

For the diffuse illumination setup, we kept the imager mounted in the same orientation but removed the set of four halogen lights used for direct illumination. Instead, we positioned two larger softbox diffuse halogen lamps (Westcott uLite) perpendicular to the snow sample surface (Fig. 3b). Aiming the lights at a large white panel on the opposite side of the snow sample and 51 cm away mimicked diffuse hemispheric illumination conditions with no direct component. Once again, we used a Spectralon panel to convert raw data to reflectance images.

Figure 4 Visible photographs of two different snow samples (a, b) are juxtaposed with greyscale images of hyperspectral reflectance at 1030 nm (c, d) as well as lidar-derived reflectance at 1064 nm (e, f). Samples are shown under direct illumination at Θ = 0° and at their original resolutions, prior to any coarsening. The data in the top row depict Sample 3, a snow microstructure with a relatively high specific surface area (26.31 mm⁻¹) compared to Sample 31 in the bottom row, which had a specific surface area of 12.40 mm⁻¹. The insensitivity of reflectance to snow microstructure in the visible range, as opposed to the dramatically different reflected magnitudes in the NIR, is apparent. The inset figures illustrate the location of the bidirectional reflectance measurements relative to the ice absorption spectra. These example spectra were produced using the Asymptotic Radiative Transfer model (Kokhanovsky and Zege, 2004).

Hyperspectral data processing

Initial processing took place in Resonon's proprietary Spectronon software, and analyses were performed thereafter in Rstudio. We began by truncating the reflectance data from each dataset to a central region of interest (ROI) encapsulating the micro-CT ROI (Fig. 3d). Our goal was to exclude edge effects from mixed pixels/points, particularly considering the lidar beam dilation when samples were tilted off-nadir. The resulting NIR-HSI ROIs contained 224 000 pixels with a spatial resolution of 0.5 mm. Reflectance images were produced from 188 of the 336 available bands, ranging from 951–1403 nm, which was trimmed to reduce noise at the lower end of the imager spectral range and to exclude regions where snow and ice are scarcely reflective. Example snow sample reflectance maps from a single band centered at 1030 nm, the location of a prominent ice absorption feature (hereafter R₁₀₃₀), are illustrated in Fig. 4c and d. The last factor we sought to examine was spatial resolution, considering that if an NIR texture signature specific to surface hoar does exist, then it is likely resolution dependent. Thus, we coarsened all datasets to spatial resolutions of 1, 2.5, 5, 10, 25, and 50 mm (hence 2 orders of magnitude) as an attempt to mimic the spatial resolutions achievable by UAV-mounted systems.

2.2.2 Lidar

Data collection

We used a Riegl VZ-6000 terrestrial laser scanner mounted to a tripod. The scanner achieves vertical (line) scanning via an oscillating mirror while moving horizontally on a rotating head (Fig. 3c). The maximum vertical scan field of view (sFOV) is 60–120° from zenith and selectable therein, while the horizontal sFOV can range from 0° to a full 360° panorama. The laser operates in the NIR range narrowly around a central wavelength of 1064 nm. We set the vertical and horizontal angular increments to 0.01° to maximize resolution (point density) and selected a pulse repetition frequency of 300 kHz. The initial laser beam diameter upon exiting the scanner is 15 mm, with a divergence of 0.12 mrad. We positioned the scanner about 2.5 m from the snow samples. The duration of each scan was approximately 1 min. The resulting data product for lidar is a “cloud” of discrete vector data points, called returns. Similar to the hyperspectral imaging setup, we used a Lambertian reflectance standard to convert the return power to the bidirectional reflectance factor for each individual point in the clouds.

Figure 5 A map of R₁₀₃₀ from NIR-HSI at its raw resolution of 0.5 mm (a) is coarsened by 2 orders of magnitude to 5 mm (b) and 50 mm (c). For each resolution, the localized standard deviation, σ₁₀₃₀, is calculated using a 9-pixel moving-window analysis (d–f). The data shown are from Sample 26 under direct illumination at an incidence angle of Θ = 0°.

Lidar data processing

Initial processing of point clouds took place in Riegl's RiScan application, and analyses were performed thereafter in Rstudio. As with the NIR-HSI data, we trimmed the point clouds to a central ROI (Fig. 3d). After truncation, the average point cloud contained 80 000 returns, with a spacing of ∼ 1.4 pts mm⁻². In order to perform pixelwise operations and to better compare with hyperspectral imagery, we converted the point cloud into an image of continuous pixels. Using a Delaunay triangulation, we interpolated point clouds into images with 1 mm spatial resolution, resulting in 56 000 pixels per dataset. The Riegl VZ-6000 lidar unit operates narrowly around a central wavelength of 1064 nm. Therefore, unlike NIR-HSI, the lidar is only capable of measuring bidirectional reflectance specifically at 1064 nm (hereafter R₁₀₆₄; e.g., Fig. 4e and f). This wavelength occurs on the shoulder of the ice absorption feature centered at 1030 nm (inset of Fig. 4e). Last, as with NIR-HSI, we coarsened all lidar datasets to the resolutions listed in the “Hyperspectral data processing” section to examine the influence of spatial resolution on NIR texture. However, we remind the reader that the original lidar beam diameter is 15 mm, and thus spatial resolutions lower than this are not as relevant as with NIR-HSI.

2.3 Texture analysis

2.3.1 Moving-window focal analysis

To restate our hypothesis, we anticipated that surface hoar would display heightened NIR variability, or texture, relative to other snow surface grain shapes due to the physical phenomenon illustrated in Fig. 1. Therefore, we sought to evaluate the localized variability of reflectance. To achieve this, we performed a moving-window focal analysis to create maps of local standard deviation. Using R₁₀₃₀ from NIR-HSI, a demonstration of both coarsening and subsequent moving-window analysis is presented in Fig. 5 for Sample 20. Beginning with a map of R at either the native (0.5 mm) resolution (Fig. 5a) or coarsened resolutions of 5 and 50 mm (Fig. 5b and c, respectively), a 9-pixel neighborhood is placed around a central pixel. The standard deviation of R₁₀₃₀ (hereafter σ₁₀₃₀) is calculated within the window via Eq. (2):

$$\begin{array}{}\text{(2)}& \mathit{\sigma }=\sqrt{{\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{3}}\sum _{j=\mathrm{1}}^{\mathrm{3}}{\left(R_{ij}-\stackrel{\mathrm{‾}}{R}\right)}^{\mathrm{2}}}{N}}},\end{array}$$

where R_ij is the reflectance value at row i and column j in the 3 × 3 grid, $\stackrel{\mathrm{‾}}{R}$, is the mean of the reflectance values in the 3 × 3 grid, and N is the total number of pixel values (nine in this case). The resulting value of σ is assigned to the central pixel.

We handled edges by truncating the window size when necessary, such that corner pixels only considered three neighboring cells, for example. Moving the window across each R₁₀₃₀ map and evaluating every pixel independently yields a map of σ₁₀₃₀ (Fig. 5d–f). Hence, these maps depict NIR reflectance texture rather than magnitude. We performed this process for all datasets – thus, for each of the 188 NIR-HSI bands and the lidar-derived R₁₀₆₄ for each illumination condition and incidence angle.

Figure 6 Flowchart of the samplewise data processing and statistical analysis workflow.

2.3.2 Spatial and spectral texture analysis

Once texture maps were derived for each dataset, we determined which spatial resolutions were best for surface hoar delineation, assuming some degree of resolution dependence. Furthermore, leveraging the spectral data provided by NIR-HSI, we examined the relationship between NIR texture and wavelength. We began by grouping the pixelwise σ distributions from surface hoar Samples 24–26 for each dataset. Similarly, we grouped the σ distributions of all other samples (hereafter termed “Other”, meaning a microstructure other than SH). Next, we calculated the median values of both grouped distributions for each instrument, band, and spatial resolution and determined the percent difference in medians for each case. Using heat maps and line graphs, we evaluated the σ distributions and the resulting differences in medians (i.e., ΔM(σ)). This allowed us to determine the spatial resolutions and, in the case of NIR-HSI, the wavelengths that maximized the difference between the SH and Other microstructure medians. We note that a larger difference in texture medians between SH and Other samples should allow for greater ease of SH classification.

2.3.3 Samplewise analysis and significance testing

We next examined how texture varied between individual samples, rather than considering SH samples against an aggregation of other microstructures (Sect. 2.3.2). As a case study, we selected a single band and spatial resolution for each instrument. For NIR-HSI, we chose the band centered at 1324 nm. The results of our spatial and spectral analysis (described in Sect. 2.3.2; results to be discussed in Sect. 3.1) indicated that this was an optimal wavelength under both direct and diffuse illumination. This is a sensible finding considering that 1324 nm sits amid a prominent ice absorption feature where reflectance is particularly sensitive to the path length of ice (i.e., the optical grain size). The lidar unit has only one band, centered at 1064 nm. We elected to optimize spatial resolutions based on the results in Sect. 3.1 as well, and thus we selected 10 mm for NIR-HSI and 5 mm for lidar. To determine if surface hoar textures are significantly larger than those of differing microstructures, we compared the distributions and median values of σ across samples. As in Sect. 2.3.2, we grouped the pixelwise σ distributions of surface hoar (Samples 24–26) for each dataset. We then compared the median value of this grouped SH distribution against samplewise medians. Specifically, we performed one-sided, one-sample t tests to assess the grouped median against the median σ of each sample. The flow chart in Fig. 6 describes the processing workflow from raw reflectance images through statistical analyses for the case study. Our hypotheses can be summarized as follows, where X describes an individual sample number:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{H}_{\mathrm{0}}:M\left({\mathit{\sigma }}^{S\mathrm{24}:S\mathrm{26}}\right)=M\left({\mathit{\sigma }}^{S\left(X\right)}\right)\\ {\displaystyle }& {\displaystyle }{H}_A:M\left({\mathit{\sigma }}^{S\mathrm{24}:S\mathrm{26}}\right)>M\left({\mathit{\sigma }}^{S\left(X\right)}\right);\mathit{\alpha }=\mathrm{0.05}.\end{array}$$

Figure 7 Heat maps depicting the percent difference in median values of σ for surface hoar samples versus all other microstructures (SH minus Other) across a variety of spatial resolutions and spectral bands (a, b). Data were acquired via NIR-HSI under diffuse (a) and direct (b) illumination conditions. In (c), a line graph illustrates the same percent difference for the lidar-derived σ₁₀₆₄, with points colored according to the same scale when an increase is observed. Vertical black lines in (a) and (b) are located at the lidar wavelength of 1064 nm for reference, while the grey vertical line in (c) is centered at zero.

2.4 Classification algorithm

Continuing our samplewise case study (Sect. 2.3.3), we next produced classified maps of surface hoar. We determined optimal threshold values of σ to delineate surface hoar from other snow surface shapes on a per-pixel basis. In addition to visualizing distributions (across all pixels) of σ with box plots for each sample, we also constructed probability density functions (PDFs). As in Sect. 2.3.2, we grouped the distributions of surface hoar samples (Samples 24–26) and compared them to the aggregated distributions of all other samples. We selected values of σ corresponding to the intersection of the grouped PDFs as the optimal thresholds of delineation for each dataset (i.e., each combination of instrument, illumination condition, and incidence angle), termed σ_crit. Using these threshold values, we performed a binary pixelwise classification: pixels with σ values above σ_crit were classified as surface hoar, while values beneath were designated as Other microstructures. We ran the classification algorithm on all samples using the appropriate σ_crit value for each dataset. To evaluate the success of the classification algorithm, we calculated the true-positive rate (TPR), true-negative rate (TNR), and overall accuracy (A) for each sample using the following equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{TPR}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{TNR}={\displaystyle \frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}}\text{(5)}& {\displaystyle }& {\displaystyle }A={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}},\end{array}$$

where TP is a true positive, TN is a true negative, etc. In this case, a TP refers to the correct identification of surface hoar when it is present, a TN corresponds to the correct identification of an Other microstructure, an FP is when the algorithm makes an incorrect surface hoar classification, and an FN is when surface hoar is misclassified as Other.

2.5 Method repeatability assessment

Our final investigation was to test the spatial mapping capacity and repeatability of our texture-based classification on a new snow sample, one that was not involved in the initial analysis. Using the techniques outlined in Sect. 2.1.1, we prepared a snow sample of rounded grains and grew surface hoar atop roughly half of the surface area, while the other half remained covered. Thus, the resulting snow surface grain shape was a 1:1 ratio of SH:RG. We proceeded to run the classification algorithm (Sect. 2.4) on the mixed sample to produce a map of surface hoar extent, using the appropriate values of σ_crit for each dataset. Unfortunately, the lidar unit was no longer available at the time of this assessment, so only NIR-HSI was evaluated. We again quantified accuracy using Eqs. (3)–(5).

3 Results

Here, we discuss results at nadir orientations, hence Θ = 0° (as defined in Fig. 3). The results of each case at Θ = 10° were very similar to their nadir counterparts, and so, for clarity, we address only the latter here. Results for Θ = 10° are presented in Appendix A.

3.1 Spatial and spectral texture analysis

We evaluated a grouped distribution of σ for surface hoar (Samples 24–26) against a grouped distribution containing all other samples for each instrument–illumination condition. We performed this evaluation across a variety of spatial resolutions spanning 2 orders of magnitude. Further, in the case of NIR-HSI where spectral data were available, we examined results over 188 bands from 951–1403 nm. Using NIR-HSI under both direct and diffuse illumination, we found that the median value of σ for surface hoar was nearly always greater than that of the Other microstructure group, but with considerable spatial and spectral dependence (Fig. 7a and b). When evaluating the difference in medians (i.e., SH minus Other), a normal distribution is evident along the spatial (vertical) axis, with the largest differences observed at 10 mm spatial resolution for both illumination cases. Spectrally, the difference in medians remains fairly constant, peaking in the ∼ 1250–1350 nm range. The maximum increase was 383 % at 1246 nm for diffuse illumination and 294 % at 1369 nm for direct illumination. The differences in median values tend to be larger in the case of diffuse illumination relative to direct, and therefore we anticipate greater ease of surface hoar classification with diffuse lighting.

To reiterate, when using NIR-HSI, median values of σ were larger in grouped SH distributions than in other microstructure distributions across nearly all cases (Fig. 7a and b); at worst, the medians were essentially equal. This provides evidence for our hypothesis that surface hoar will produce increased NIR reflectance texture under a variety of conditions. Our lidar observations, however, were less consistent. The mere presence of heightened texture in SH was dependent on spatial resolution (Fig. 7c). At the lowest (1 and 2.5 mm) and highest (50 mm) spatial resolutions, we observed larger medians in the Other microstructure group than in SH, a result that contrasts with our hypothesis. A normal distribution of delta median values corresponding to spatial resolution is evident, as in NIR-HSI, with the 5 and 10 mm datasets producing the largest (positive) differences. As discussed, we were limited to the evaluation of a single band with our lidar unit (R₁₀₆₄), although our NIR-HSI spectral results (Fig. 7a and b) indicate that 1064 nm is a suitable wavelength.

Figure 8 Samplewise median values of reflectance (a–c) and σ (d–f) for each instrument–illumination case study. The black lines are spline or linear best fits to all samples other than surface hoar, while fuchsia horizontal reference lines depict SH median values. Colors and point shapes correspond to those suggested by the ICSSG (Fierz et al., 2009).

3.2 Samplewise statistical analysis and significance testing

We conducted a samplewise case study using the NIR-HSI-derived σ₁₃₂₄ (both direct and diffuse) at 10 mm and the lidar σ₁₀₆₄ at 5 mm spatial resolution. These selections were based on our findings in Sect. 3.1; the spatial resolutions were the optimal choice for each instrument, and the 1324 nm maximized NIR-HSI texture increases under both illumination conditions. When using NIR-HSI, we found that surface hoar exhibited larger values of σ₁₃₂₄ relative to other sample microstructures in both direct and diffuse illumination conditions. This is evident in Fig. 8d and e. Figure 8 also outlines the problem with using NIR reflectance magnitude to delineate surface hoar. In both cases, the median reflectance magnitude (R₁₃₂₄) gradually increased proceeding from lower to higher SSA (Fig. 8a and b), as expected. This leaves the reflectance magnitude of surface hoar “hidden” in the middle of all other grain shapes, with moderate median R₁₃₂₄ values (horizontal fuschia reference lines). When examining values of σ₁₃₂₄ (Fig. 8d and e), a different pattern is evident. In Fig. 8d and e, median values of σ₁₃₂₄ are fairly constant across all samples with the exception of surface hoar, where a spike is present. The separation between the median value of σ₁₃₂₄ for surface hoar samples compared to the best-fit line for all other samples (black lines) further illustrates this and may allow for the delineation of surface hoar based on a texture parameter. Our results for lidar R₁₀₆₄ indicate a similar trend, although the distinction is less pronounced (Fig. 8c and f).

Figure 9 Samplewise box-plot analysis for all three instrument–illumination scenarios. Boxes are colored by primary grain shape following the ICSSG. Horizontal fuschia reference lines depict the median value of σ for surface hoar Samples 24–26. Samples with an asterisk had median values significantly lower than surface hoar medians based on one-sample, one-sided t tests at α=0.05.

Table 2 Sample median values of σ for each instrument and illumination condition. Median σ values of surface hoar Samples 24–26 are shown at the top. Bold sample median σ values are significantly lower than the corresponding surface hoar medians at a significance level of α=0.05.

To restate our central finding, when using NIR-HSI under both direct and diffuse illumination, surface hoar exhibited larger median values of σ₁₃₂₄ relative to other sample microstructures in all cases (Table 2, Fig. 9a and b), thus confirming our hypothesis. Furthermore, this increase was statistically significant in all cases. The in situ Samples 25 and 26 were particularly distinct, with interquartile ranges rarely overlapping those of any other sample, making these easily discernible. Interestingly, Sample 24, which consisted of sieved SH grains (Table 1), displayed heightened texture under direct but not diffuse illumination. As with Fig. 8, our lidar results were similar to those of NIR-HSI but less pronounced (Fig. 9c), with particular confusion occurring with PP and MF samples. Still, the grouped SH median value of σ₁₀₆₄ was significantly larger than the sample median in $\mathrm{30}/\mathrm{38}$ cases. Table 2 provides complete significance findings for all scenarios, with greyed cells denoting a case where the sample median σ value is significantly less than the surface hoar median.

Figure 10 Probability density functions juxtaposing the σ distributions of surface hoar Samples 24–26 (fuchsia curves) with those of all other samples (grey curves). The area under each curve is equal to unity. Dotted vertical reference lines represent the distribution intersection, where the critical-threshold values (σ_crit) were extracted. Median accuracy metrics from the resulting samplewise binary classifications are also listed for each scenario.

Figure 11 Example transformations from maps of σ₁₃₂₄ (a–f) to pixelwise classifications (g–l) via binarization with a critical threshold. The data shown here are from direct illumination at 10 mm spatial resolution.

Table 3 Samplewise and median accuracy (A) values for each instrument–illumination combination. Results for surface hoar Samples 24–26 are italicized.

3.3 Classification algorithm

We used PDFs to compare the probability density distributions of surface hoar samples (Samples 24–26) against all other samples, which allowed us to determine critical texture thresholds (σ_crit) for surface hoar delineation under each instrument–illumination scenario. This process is illustrated in Fig. 10. Probability density describes the relative likelihood of sampling a certain value of σ; y values are assigned such that the entire area under the curve is equal to unity. Therefore, the area under the curve over a given range of σ is equivalent to the probability of sampling a value of σ in that range. Selecting the intersection point of the two distributions, following Champollion et al. (2013), allows for optimal binary classification of a given pixel. The resulting critical values for each condition are listed on the PDF plots. Using these σ_crit values, we conducted the pixelwise classification of all samples). Examples of this binarization for several samples of varying primary grain shape are shown in Fig. 11. All samplewise accuracy values (A) are listed in Table 3, while median values of A, TPR, and TNR for each scenario are included on the PDF plots in Fig. 10. Classification results generally mirrored those of the statistical analysis. That is, results were excellent for NIR-HSI under both direct and diffuse illumination but dwindled substantially when using the lidar-derived σ₁₀₆₄. The classification algorithm generally struggled most with samples with a FC or MF primary grain shape.

Figure 12 A sample with a 50:50 ratio of RG : SH surface grain shapes is displayed in the visible (a). Under both diffuse (b, d, f) and direct (c, e, g) illumination conditions, the bidirectional reflectance factor at 1324 nm is extracted (b, c) from NIR-HSI. Via a moving-window analysis, localized standard deviation (σ₁₃₂₄), or NIR texture, is quantified (d, e) and then binarized using critical thresholds to produce classified maps (f, g). Overall accuracy, true-positive rates, and true-negative rates are listed on the classified maps. The transitional zone was excluded in accuracy analyses.

3.4 Classification assessment

Classification accuracy, assessed on a new sample comprising a 1:1 ratio of RG and SH surface grain shapes, proved consistent, demonstrating repeatability of the texture phenomenon when using NIR-HSI (Fig. 12). Further, we observe the capacity to map surface hoar extent amid mixed surface conditions. In the R₁₃₂₄ maps (Fig. 12b and c), the heterogeneity of surface hoar reflectance is already apparent compared to the RGs, and this was quantified via a moving-window analysis (Fig. 12d and e). Using the appropriate values of σ_crit (Sect. 3.3), we binarized each texture image to create a classified data product (Fig. 12f and g). In the classified maps, the rate of correct rejection of the RG surface (i.e., the TNR) was perfect (100.00 %) for diffuse illumination and suitable (92.23 %) under direct illumination. The rate of accurate identification of SH (i.e., the TPR) was also excellent, at 99.35 % and 99.32 % for diffuse and direct illumination, respectively, while the overall accuracy was 99.61 % and 96.15 %.

4 Discussion

4.1 NIR-HSI

We found that, when measured with NIR-HSI, surface hoar exhibited larger median values of σ, an NIR texture metric, than other snow microstructures across a variety of spatial resolutions and spectral bands. Furthermore, in a samplewise case study evaluating σ₁₃₄₂ at 10 mm spatial resolution, we determined that median surface hoar values were significantly larger than other snow surface structures (the primary grain habits of PP, DF, FC, RG, and MF) in all cases under both direct and diffuse illumination (one-sided, one-sample t test, α=0.05). Our findings are consistent with Champollion et al. (2013), who found that under artificial lighting conditions, the NIR texture (in this case, a contrast index) of the Antarctic snow surface was higher when surface hoar was present as compared to when it was absent. These researchers used a passive NIR camera and downward-looking lights, all mounted ∼ 2 m above a flat snow surface, to create a field setup akin to the direct-illumination (Θ = 0°) case presented here. Although they do not explicitly mention the spatial resolution of their images, they imply that individual surface hoar crystals spanned 5–10 pixels, and thus their resolution was likely on the order of millimeters. While Champollion et al. (2013) simply juxtaposed the presence versus absence of surface hoar, we build on those foundational findings by quantifying the texture phenomenon against a variety of differing, thoroughly characterized snow microstructures in a controlled laboratory environment. Further, our work features the explicit inclusion of diffuse illumination, varied illumination/viewing incidence angle, and a range of spatial resolutions and spectral bands. Last, whereas Champollion et al. (2013) used texture metrics across an entire image to confirm the presence of hoar crystals on a given day (with an accuracy of 94 %), we extend this classification by demonstrating a per-pixel mapping methodology, allowing us to map the spatial extent of surface hoar within an image (Fig. 12). Median sample classification accuracy was 96.91 % and 97.37 % for diffuse and direct illumination, respectively. As discussed in Sect. 3.1 and illustrated in Fig. 8, surface hoar exhibits moderate NIR reflectance when compared to other grain shapes, and thus leveraging reflectance magnitude alone is insufficient for surface hoar delineation. These results highlight the importance of our texture-based approach.

Spatial resolution was perhaps the factor with the most influence on our results, with the maximal difference between surface hoar and other median values occurring at 10 mm resolution for both illumination cases (Fig. 7). However, consistent differences in median values can be observed at nearly all resolutions beneath 50 mm, indicating that classification at coarser resolutions may still be suitable under the right conditions. This trend makes sense physically; we attribute the rise in texture metrics associated with surface hoar as being related to the lateral spacing between large hoar crystals (Fig. 1). Given that this spacing is often on the order of millimeters, it is possible that the raw resolution is too fine to optimally observe this variability, while 50 mm resolution is too coarse. Further, the fact that results were fairly consistent spectrally from 951–1403 nm implies that hyperspectral capacity is likely unnecessary; a broadband passive NIR sensor could likely observe the same texture increases. An even more cost-effective option would be to remove the NIR filter from a standard CCD camera, allowing for light detection up to 950 nm. The success of Walter et al. (2023), discussed further in Sect. 4.2, provides optimism for observing SH texture at NIR wavelengths beneath 950 nm. The samples that proved the most challenging to discern from SH were FC (Samples 19–23) and the large MF grains from Batch L (Samples 37–41). The latter was likely due to enhanced surface roughness, as these grains were several mm in length (Fig. 2) – a rather extreme case, and this indicates that caution should be used when conditions favor wet snow metamorphism. Sporadic misclassification of FC (e.g., Fig. 11i) is perhaps understandable, given that these grains form from kinetic snow metamorphism, which is somewhat similar to surface hoar growth. In practice, it may be useful to try to identify these surfaces as well, as near-surface FC also tend to act as weak layers once buried. The relatively low median σ₁₃₄₂ value of the sieved SH Sample 24 under diffuse illumination (Fig. 9a) is an interesting result. This anomaly perhaps implies that the texture signature dissipates when the predominately vertical orientation of SH grains is interrupted, and thus surface hoar that has been blown over is likely more difficult to detect. The low accuracy value of this sample under diffuse illumination (47.76 %, Table 3) is evidence of this. The uncertainty involved in SH classification for any case can be quantified by examining the area under the intersecting PDF curves (Fig. 10). Furthermore, uncertainty with regards to a specific microstructure (such as FC) could be observed by comparing PDFs between individual samples, rather than by grouping SH samples versus all other samples. However, our goal was to determine robust thresholds for the delineation of SH from any other microstructure.

4.2 Lidar

As impressive as the performance of NIR-HSI was, the results of our lidar analysis were more perplexing. In Fig. 7c, we can observe that the median value of σ₁₀₆₄ for SH Samples 24–26 was only greater than that of other microstructures at spatial resolutions of 5, 10, and 25 mm. At the largest and smallest resolutions, much like with NIR-HSI, performance diminished. However, unlike NIR-HSI, in these cases, the median SH value was actually lower than the other microstructure median, a result that runs counter to our hypothesis. Even at 5 mm resolution, our lidar classification results were considerably less robust than those of NIR-HSI, with a median sample accuracy of 66.99 %. Unexpected results in texture analyses are not unheard of. While Champollion et al. (2013) noted increases in NIR texture when surface hoar was present under direct “artificial” illumination at nadir incidence (as we did with NIR-HSI), they also documented a reversal of this trend under natural (solar) illumination. However, this was at very large (80–85°) solar zenith angles and thus represents an extreme case of grazing incidence, and the cause of this observation was not thoroughly explained. For the laboratory setup used here, we expected lidar to reproduce similar results to those of our NIR-HSI analysis and the findings of Champollion et al. (2013) using artificial light. Like NIR-HSI, the lidar classification algorithm struggled with the large MFs of Batch L as well as PP (Fig. 9c).

In general, the use of lidar reflectance to ascertain snow surface properties is far less validated relative to passive imagery, like NIR-HSI. While passive NIR imagery has been used to estimate SSA or r_e for decades, leveraging the well-established sensitivity of NIR reflectance to snow microstructure, very few studies (e.g., Yang et al., 2017) have attempted to do so with lidar. Therefore, the dependence of NIR lidar reflectance on snow microstructure is substantially less understood. This is important because lidar scanning represents a unique bidirectional reflectance scenario. For instance, the beam emitted from lidar units is collimated, and thus it experiences less loss due to scattering and absorption compared to direct solar irradiance and results in higher irradiance at the wavelength of interest. Collimation alters scattering and reflectance in optically rough materials like snow relative to an un-collimated illumination source (Murphy, 2006). Additionally, lidar beams are predominately linearly polarized, which contrasts with the often randomly polarized nature of solar illumination (Sassen, 2005). The scattering process is partly polarizing, and therefore multiple scattering can have a significant polarization dependence (Li et al., 2008; Bhandari et al., 2011). Last, lidar presents a monostatic geometric condition on bidirectional reflectance such that every measurement strictly observes direct backscatter. This is beneficial in that it limits the number of characterizations required compared to a full bidirectional reflectance investigation. However, observations of snow reflectance at this geometry are lacking, and radiative transfer models are not well-validated for the case of direct backscatter. Thus, lidar is an optically unique case of bidirectional reflectance that requires careful examination, so it is not entirely surprising that a texture analysis of lidar reflectance produced peculiar results.

One possible explanation for the relatively poor performance of our algorithm when using lidar reflectance can be found in the work of Walter et al. (2023), who leveraged a laboratory setup much like the apparatus presented here. These researchers focused on a temporal analysis, using a 905 nm lidar to observe changes in the mean reflectance magnitude and standard deviation (across an entire image) at prescribed time intervals as they grew surface hoar atop a layer of compressed PP. Consistent with our spatial analysis, they noted that reflectance magnitude was insufficient to characterize surface hoar growth, as the mean reflectance changed by only 4 %, while standard deviation increased by as much as 600 %. Although this juxtaposition features only a single microstructure other than SH (compressed PP), these lidar-based results seem more encouraging than our own and are more consistent with our NIR-HSI findings. A key distinction is, perhaps, the lidar spot size. Their spot size of 3 mm is much smaller than that of our lidar setup (15 mm) and likely closer to the size and spacing of individual surface hoar crystals in many cases. Indeed, the authors predict that surface hoar crystals smaller than the lidar beam diameter will not be detectable, and thus it is possible that our beam diameter is simply too coarse for this application.

4.3 Future work and speculations on scaling to field applications

Though scaling to field applications presents considerable challenges and uncertainty, it is likely that our NIR-HSI findings can be extended to operational avalanche forecasting in the near future. A simple setup like that of Champollion et al. (2013), where a downward-looking NIR camera acquired daily and nightly images, could be installed at remote weather stations to monitor the formation and persistence of surface hoar layers prior to burial. Such remote measurements are not currently available but are critically important for avalanche forecasters. Further, UAV snow mapping using NIR-HSI has recently been demonstrated as a useful tool to measure snow grain size and albedo at the slope scale (Skiles et al., 2023), and future studies should consider texture analyses. To accomplish this, a greater variety of incidence angles must be evaluated. Although it is encouraging that our results remained promising at Θ = 10° (Appendix A), more oblique angles will inevitably be encountered in the field. Even on overcast days, when diffuse solar illumination is prevalent, variance in the terrain slope beneath a downward-looking imager would still alter the viewing incidence angle and thus the magnitude, and likely texture, of NIR snow reflectance. Using direct solar illumination would add another factor, the illumination incidence angle, although this could be kept consistent by using artificial lighting. Ideal thresholds of σ_crit will likely continue to vary between incidence angles as well as between instruments/applications. Further, the flight plan and/or optical logistics would need to be controlled to ensure adequate spatial resolution. Another factor that may be worth investigating is the window size during focal analysis, although a preliminary study determined that the neighborhood size was of little consequence for our data. Last, evaluating the influence of the underlying substrate grains on SH texture would be useful, as we only examined rounded grains here.

Although we realized limited success with lidar, the use of lidar for surface hoar mapping via texture requires further evaluation, particularly with regards to beam diameter. A better physical understanding regarding the resolution dependence of our lidar results is needed. At a minimum, lidar could be useful in conjunction with NIR-HSI or another passive NIR detection method. Using the spatial capacity of lidar, values of slope angle (and thus incidence angle if using a downward-oriented imager) could be determined on a per-pixel basis, allowing for the proper threshold value to be used when analyzing NIR-HSI texture data. Though challenging, slope- or even basin-scale mapping of surface hoar extent is likely possible with current technology, though a field campaign would be required to tune our method to a wider variety of illumination conditions and incidence angles. Maps of surface hoar extent over avalanche terrain prior to burial would be vital in making slope-specific avalanche forecasts. Such maps could also identify likely avalanche trigger points, improving avalanche mitigation with explosives, for example. While our work focused on surface hoar classification, texture analysis will likely provide a useful tool for evaluating other physical phenomena as well, such as surface roughness.

5 Conclusion

Our research demonstrates a novel method for mapping the spatial extent of surface hoar using NIR texture in a cold laboratory. In essence, we found that:

• i.

  Hyperspectral imaging can robustly measure the texture of snow and ice by computing the pixelwise variability in reflectance at any NIR wavelength.

• ii.

  When evaluated with NIR-HSI, surface hoar has significantly heightened NIR texture relative to other snow microstructures across a range of spectral bands and spatial resolutions, likely as a result of variable ice absorption and specular contributions.

• iii.

  Near-infrared texture thresholds can be used to binarize NIR-HSI texture measurements, resulting in accurate maps of surface hoar spatial extent on a per-pixel basis.

• iv.

  Similar results were achieved with 1064 nm lidar, although the phenomenon was resolution dependent and the performance was substantially less robust. The use of lidar for this purpose requires further investigation, and is likely dependent on beam diameter.

As NIR-HSI and lidar become more economical, these may provide capable methods for measuring NIR texture and subsequently mapping surface hoar. Extending the work presented here to field operations will have immediate implications for broad-scale snow surface mapping and avalanche forecasting.

Appendix A

Results for the Θ = 10° incidence viewing angle case are presented below. As mentioned in the text, the results are very similar to the case of Θ = 0°, particularly regarding the samplewise comparison and classification mapping. This is perhaps unsurprising, as the bidirectional reflectance distribution function (BRDF) of snow changes very little for the case of direct backscatter between 0–10°. Thus, we certainly would not expect large changes in reflectance magnitude, and this is likely true for spatial variability as well. We observed the most pronounced differences in the spatial/spectral analysis. By juxtaposing Figs. 7 and A1, we note slightly larger increases in M(σ) under diffuse illumination for Θ = 10° relative to Θ = 0°, while the opposite is true under direct illumination. The latter difference is likely related to how non-nadir direct illumination interacts with the predominately vertically oriented SH crystals, although this topic requires further investigation. As mentioned in Sect. 4.3, it is imperative that future studies evaluate texture at more oblique incidence angles, particularly when considering scaling to field applications.

Figure A1 Heat maps depicting the percent difference in median values of σ for surface hoar samples versus all other microstructures (SH minus Other) across a variety of spatial resolutions and spectral bands (a, b). Data were acquired via NIR-HSI under diffuse (a) and direct (b) illumination conditions. In (c), a line graph illustrates the same percent difference for the lidar-derived σ₁₀₆₄, with points colored according to the same scale when an increase is observed. Vertical black lines in (a) and (b) are located at the lidar wavelength of 1064 nm for reference, while the vertical grey line in (c) is centered at zero. Relative to the nadir case, performance generally improved for NIR-HSI under diffuse illumination but decreased under direct illumination for both NIR-HSI and lidar.

Figure A2 Samplewise median values of reflectance (a–c) and σ (d–f) for each instrument–illumination case study. The black lines are spline or linear best fits to all samples other than surface hoar, while horizontal fuchsia reference lines depict SH median values. Colors and point shapes correspond to those suggested by the ICSSG (Fierz et al., 2009).

Figure A3 Samplewise box-plot analysis for all three instrument–illumination scenarios. Boxes are colored by primary grain shape following the ICSSG. Horizontal fuchsia reference lines depict the median value of σ for surface hoar Samples 24–26. Samples with an asterisk had median values significantly lower than surface hoar medians based on one-sample, one-sided t tests at α=0.05.

Table A1 Sample median values of σ for each instrument and illumination condition. Median σ values of surface hoar Samples 24–26 are shown at the top. Bold sample median σ values are significantly lower than the corresponding surface hoar medians at a significance level of α=0.05.

Figure A4 Probability density functions juxtaposing the σ distributions of surface hoar Samples 24–26 (fuchsia curves) with those of all other samples (grey curves). Dotted vertical reference lines represent the distribution intersection, where the critical-threshold values (σ_crit) were extracted. Median accuracy metrics from the resulting samplewise binary classifications are also listed for each scenario.

Table A2 Samplewise and median accuracy (A) values for each instrument–illumination combination. Results for surface hoar Samples 24–26 are italicized.