Mapping liquid water content in snow: An intercomparison of mixed-phase optical property models using hyperspectral imaging and in situ measurements

It is well understood that the distribution and quantity of liquid water in snow is relevant for snow hydrology and avalanche forecasting, yet detecting and quantifying liquid water in snow remains a challenge from the microto the macro-scale. Using near-infrared (NIR) spectral reflectance measurements, previous case studies 10 have demonstrated the capability to retrieve surface liquid water content (LWC) of wet snow by leveraging shifts in the complex refractive index between ice and water. However, different models to represent mixed-phase optical properties have been proposed, including (1) internally mixed ice and water spheres, (2) internally mixed water coated ice spheres, and (3) externally mixed interstitial ice and water spheres. Here, from within a controlled laboratory environment, we determine the optimal mixed-phase optical property model for simulating wet snow 15 reflectance using a combination of NIR hyperspectral imaging, radiative transfer simulations (DISORT), and an independent dielectric LWC measurement (SLF Snow Sensor). Maps of LWC were produced by finding the least residual between measured reflectance and simulated reflectance in spectral libraries, generated for each model with varying LWC and grain size, and assessed against the in situ LWC sensor. Our results show that the externally mixed model performed the best, retrieving LWC with an uncertainty of ~1%, while the simultaneously retrieved 20 grain size better represented wet snow relative to the established scaled band area method. Furthermore, the LWC retrieval method was demonstrated in the field, imaging a snowpit sidewall during melt conditions, mapping pooling water, flow features, and LWC distribution in unprecedented detail.


Introduction
The distribution and quantity of liquid water within a snowpack, introduced by rain and/or melt, is relevant for 25 multiple snow related applications including snow hydrology, remote sensing, and avalanche forecasting. In terms of snow hydrology, water is an indicator of snow energy balance and snow melt timing as the change in phase from ice to water indicates that the cold content of the snowpack is depleted, and that energy balance inputs are contributing to melt (Dewalle and Rango, 2008). Rain-on-snow can accelerate this process by contributing large energy inputs into the snowpack over a short amount of time (Mazurkiewicz et al., 2008). Water at the surface will also lower 30 snow albedo, initiating a positive feedback loop that increases absorbed solar radiation, the main driver of snowmelt (Gupta et al., 2005). For active and passive microwave remote sensing of snow, the presence of water alters microwave signatures because of the large difference in relative permittivity between liquid water and ice (i.e., dry snow). For active microwave sensors, wet snow causes characteristic changes in microwave backscatter and reduces https://doi.org/10.5194/tc-2021-247 Preprint. Discussion started: 12 August 2021 c Author(s) 2021. CC BY 4.0 License. penetration depth (Shi and Dozier, 1992), while for passive sensors, the emissivity of the snow surface is increased 35 (Walker and Goodison, 1993). For avalanche forecasting, the infiltration of liquid water into the snowpack impacts snow stability (Conway and Raymond, 1993). The strength of the snowpack can be increased at lower water content, where grains form well bonded clusters, but reduced at higher water content when water flow through pore space deteriorates a significant number of snow grain bonds resulting in relatively cohesionless particles (Colbeck, 1982).
Although it is recognized as a critical snow property across the cryospheric sciences, liquid water content (LWC) 40 measurements in a snowpack are notoriously difficult to accurately quantify due to the high spatial and temporal variability of liquid water distribution.
Here, the utility of mapping LWC in situ using near-infrared hyperspectral imaging (NIR-HSI) and radiative transfer model inversion is assessed. This approach leverages the segments of the near-infrared (NIR) spectrum where the optical properties of liquid water, hereafter referred to as water, vary from those of ice. To date, wet snow 45 has been modeled using effective spheres with a known radius, referred to as the effective grain radius (re), where the optical properties of ice and water are mixed either internally or externally. In wet snow, the arrangement of water relative to ice particles and pore space varies based on the level of saturation, which may be relevant for radiative transfer modeling. For example, water saturation below 7% (pendular regime), when water is contained in menisci held in between the ice particles (Colbeck, 1979), might be best represented using an internally mixed 50 particle model where an ice sphere is coated in water. On the other hand, at saturation above 7% (funicular regime), when ice particles become surrounded by water as it fills the pore space, this might be best represented as an internally mixed-phase sphere or externally mixed media model using interstitial ice and water spheres. Although different mixing model representations have been proposed and demonstrated, no study has quantitively compared the different approaches, or compared LWC retrievals to established LWC measurement 55 methods. Without intercomparing or validation, the best approach for retrieving LWC from NIR spectral reflectance has yet to be determined. Additionally, radiative transfer approaches to retrieving re are based on the optical properties of ice and implicitly assume dry snow, and such retrievals have not been assessed for wet snow. The main objectives of this study are threefold: (1) intercompare three wet snow reflectance models against measured LWC from a dielectric measurement instrument in a controlled laboratory environment, (2) simultaneously assess 60 effective grain size retrieval methods and their suitability for use with wet snow, and (3) demonstrate the capability of a compact NIR hyperspectral imager to map LWC and snow grain size at the laboratory and field scales.

Liquid water in snow
Water movement through snow is a spatially and temporally complex process, controlled by water saturation level, 65 snow microstructure, and snow topography. At low water saturation (<7%), capillary forces hold water in the pore space between grains, resulting in a relatively slow and semi-uniform horizontal infiltration front, referred to as matrix flow (Colbeck, 1973). On the other hand, at high saturation (>7%), water can flow through the snowpack preferentially, resulting in the formation of vertical "finger" flows at the centimeter scale or macropore flows at the meter scale (Schneebeli, 1995;Newman et al., 2004). Although gravitational forces primarily drive vertical movement of water in snow, large amounts of water can be diverted horizontally due to stratigraphic layers in the snowpack, such as ice crusts or capillary barriers (i.e., fine grains over coarse grains) (Waldner et al., 2004;Webb et al., 2021).

Measurement of liquid water in snow
The complexity of water movement through snow makes observations and measurements challenging. Early 75 observations of water flow patterns through snow were made using dye tracers (Seligman et al., 1936;Gerdel, 1954), a method which is still used today. Dye tracers provide a spatial visualization of water infiltration that has been used to study processes such as preferential flow (Schneebeli, 1995) and capillary barriers (Avanzi et al., 2016). Additionally, fluorescent dye tracers have been used to classify LWC regimes across a snowpit wall (Waldner et al., 2004), however, the method remains primarily a qualitative visualization technique.

80
In situ measurements of LWC in snow have traditionally been measured by centrifugal separation (Kuroda, 1954), melting calorimetry (Yosida, 1940), freezing calorimetry (Jones et al., 1983), and the dilution method (Davis et al., 1993). A more detailed summary of these methods can be found in Stein et al. (1997). Generally, these methods are difficult to perform, time consuming, and have only been occasionally used since their introduction.
More commonly, LWC measurements are obtained using dielectric methods at frequencies ranging from 1 MHz to 1 85 GHz by leveraging the large differences in the relative permittivity (εr) between water (εr ≈ 88), ice (εr ≈ 3.15), and air (εr ≈ 1) . This is done by time domain reflectometry (TDR) or with capacitance sensors which measure the relative permittivity of snow (Lundberg, 1997;Denoth et al., 1984). Although the measured relative permittivity is primarily a function of LWC, snow density also has some influence. Therefore, dielectric methods also require a separate density measurement. Examples of currently available dielectric instruments include the 90 Snow Fork (Sihvola and Tiuri, 1986), Denoth Meter (Denoth, 1994, A2 Photonics WISe Sensor (A2 Photonic Sensors, 2019), and the SLF Snow Sensor (FPGA Company, 2018), which was used in this study. Although these instruments make measurements quicker relative to traditional methods, they often require destructive sampling, and only provide a discrete volume-averaged point measurement. Therefore, there is currently no in situ method to effectively quantify spatial variability of LWC at a high (sub-cm) spatial resolution.

95
Previous non-destructive measurements of LWC in snow have been made using remote sensing techniques.
Like dielectric sensors, active and passive microwave sensors leverage the difference in relative permittivity between water, ice, and air. At the ground-based scale, upward-looking ground penetrating radar (upGPR) has been used to measure the volumetric LWC directly above antennas buried below a snowpack (Schmid et al., 2014). At the spaceborne scale, active and passive microwave sensors have been used to make classification maps of wet or dry 100 snow at spatial resolutions on the order of tens of meters (e.g., (Lund et al., 2020;Walker and Goodison, 1993)).
Similarly, in the optical wavelengths, the shift in absorption patterns of ice and water across the NIR have been leveraged to map surface LWC (Green et al., 2002), which is the primary method of interest in this work.

Modeling wet snow near-infrared reflectance
Absorption in the optical wavelengths is described by the imaginary part of the complex refractive index. Although 105 the absorption patterns across the NIR are similar between ice and water, there are shifts that distinguish the https://doi.org/10.5194/tc-2021-247 Preprint. Discussion started: 12 August 2021 c Author(s) 2021. CC BY 4.0 License. different phases. The spectral complex refractive index for ice (Warren and Brandt, 2008) and water at 0 o C (Rowe et al., 2020) across NIR wavelengths is shown in Figure 1. Compared to the difference in the relative dielectric properties across the radio and microwave wavelengths, the shifts in the imaginary part of the complex refractive index are relatively minor, and therefore require measured reflectance at multiple wavelengths and radiative transfer 110 modeling to detect. Additionally, the penetration of light in the NIR wavelengths is relatively shallow, limiting detection of water to the snow surface (~2 cm). This has limited the use of optical methods to detecting surface water using either in situ spectrometer measurements, or airborne imaging spectrometer measurements (Green et al., 2002;Hyvarinen and Lammasniemi, 1987).
Inversion of radiative transfer models is commonly used in remote sensing applications to retrieve physical 115 snow properties from measured spectra and many modeling approaches have been proposed. Hyvarinen and Lammasniemi (1987) modeled the reflectance of wet snow using a collection of spheres with radius re, to simultaneously estimate LWC and grain size. To describe the optical properties of the effective spheres, an effective complex refractive index (keff) was calculated by volume-mixing the complex refractive index of ice and water.
Using a forward modeling approach, LWC and re were retrieved using only three bands (1030, 1260 and 1370 nm), 120 which were assessed using a dilatometer in a laboratory. Alternatively, Green et al. (2002) modeled wet snow using two approaches: 1) as a collection of water coated ice spheres and 2) water spheres interspersed in the interstitial space within an ice-sphere matrix. LWC and re were retrieved by matching the simulated spectra to measured spectra by finding the lowest residual. By visual inspection, Green et al. (2002) concluded that wet snow is best modeled as water coated ice spheres, though no quantitative retrieval assessment was performed. More recently, a 125 three band ratio method to classify wet or dry snow was proposed by Shekhar et al. (2019), based on a correlation between field spectrometer measurements and Snow Fork LWC measurements, although effective grain radius was not considered. To date, three radiative transfer approaches for simulating the reflectance of wet snow have been proposed: (1) mixed phase spheres, hereafter referred to as "keff spheres", (2) water coated ice spheres, hereafter referred to as 130 "coated spheres", and (3) externally mixed interstitial ice and water spheres, hereafter referred to as "interstitial spheres". A schematic of each model is presented in Figure 2. The keff and coated sphere models are referred to as internally mixed particles because the optical properties are mixed inside of a single particle having a single re. The interstitial sphere model is referred to as an external mixture because the components are assumed to be physically separate from one another.

Methodology
Three optical property mixing models were used to simulate the bidirectional reflectance of wet snow across a range of re and LWC using radiative transfer modeling. To determine the optimal mixing model for retrieving LWC from NIR reflectance, snow samples were prepared in a controlled laboratory environment with varying grain type, grain size, and density, and then subjected to warm air advection to induce melt. During melt, time series NIR-HSI 140 measurements were taken and LWC was retrieved in each pixel by best matching the measured spectra to the simulated spectra from each of the three models. For comparison to the NIR-HSI LWC retrievals, time series LWC measurements were taken with the SLF Snow Sensor, a dielectric measurement instrument, and compared to the average LWC from pixels covering the same measurement area as the SLF Snow Sensor, such that the two measurements would be comparable on the same spatial scale. Time series measurements of re were retrieved 145 simultaneously with LWC and were compared against the established scaled band area method of Nolin and Dozier (2000), which has been previously applied to the same compact hyperspectral imager used in this study to map the re of dry snow (Donahue et al., 2021). Lastly, retrievals are demonstrated in the field across an image of a snowpit wall, visualizing water infiltration and quantifying vertical LWC and re distributions.

Near-infrared hyperspectral imager
Snow reflectance in the NIR was mapped with a Resonon Inc. Pika NIR-320 near-infrared hyperspectral imager. A brief description of the instrument follows, for a more detailed description see Donahue et al. (2021). The imager has a spectral resolution of 4.9 nm, measuring 164 bands across the NIR region from 900 -1700 nm. The imager constructs a 2-D image containing the full spectrum in each pixel by collecting the image line-by-line, known 155 commonly as a "push-broom" or "line-" scanner. Thus, to collect an image, the camera needs to be moving (translating or rotating) relative to the scene, or the scene needs to be moving relative to the imager. Here, both types of image acquisition techniques are used. In the laboratory, a linear scanning stage was used to move the sample beneath the sensor, while in the field, a rotational stage mounted on top of a tripod was used to scan the snow pit face.

SLF Snow Sensor
The SLF Snow Sensor (FPGA Company, 2018), hereafter referred to as the "SLF sensor", is a capacitance sensor that is placed on the snow surface to measure the relative permittivity. This is used to determine snow density and LWC, in dry snow and wet snow conditions, respectively. The factory calibration for the LWC measurement is based on an empirical equation derived from reference measurements of snow with varying wetness and density 165 using the dilution method (Davis et al., 1985) and weighted volumes. The sensor measures a snow surface area of 45 x 95 mm and the penetration of the electric field into the snow is ~17 mm, such that the SLF sensor produces a spatially comparable measurement to the retrieved LWC from NIR-HSI method presented here because the penetration of NIR light into snow is similarly shallow. Additionally, time series LWC measurements over the same surface area can be made because the sensor is non-destructive to the snow surface.

175
This quasi-monostatic configuration, results in a bidirectional reflectance measurement in each pixel of the image when calibrated using a white reference panel. A large Spectrolon © 99% reflectance panel was placed at the same height as the surface of the snow samples and filled the imagers entire field of view, such that each snow sample could be calibrated from radiance to reflectance on a pixel-by-pixel basis. This method of calibration is ideal for hyperspectral imaging because it minimizes effects to illumination imperfections across the scene.

180
Snow samples were prepared in the laboratory using laboratory-made and collected natural snow to generate a dataset with a range of grain types including precipitation particles (PP), decomposing and fragmented precipitation particles (DF), rounded grains (RG), melt forms (MF), and faceted crystals (FC) (Fierz, 2009). The dry snow density, measured by weighing the sample container with a known volume and using the SLF sensor, ranged between 115 and 510 3 . Detailed properties for each snow sample are given in Table 1. Snow samples were made 185 by sieving snow into a rectangular wooden sample container having dimensions of 17.1 cm x 12.4 cm x 8.5 cm (H x W x D). The wooden sample container was subdivided into three regions of interest (ROI), each having dimensions slightly larger than the SLF sensor, such that LWC measurements would not be impacted by edge effects from the sample container. The surface of each snow sample was scraped with a crystal card to create a flat surface level with the top of the sample container and to minimize surface roughness. Snow samples were then kept in a cold room at -190 10 o C for 24 hours to equilibrate and ensure the sample was completely dry (i.e., no water). https://doi.org/10.5194/tc-2021-247 Preprint. Discussion started: 12 August 2021 c Author(s) 2021. CC BY 4.0 License.
For each snow sample, an initial image was taken while the cold room was at -10 o C to obtain the dry snow condition (i.e., LWC = 0%). The dry snow surface density in each ROI was measured using the SLF sensor and used as the input dry snow density to the SLF sensor for the subsequent LWC measurements. Following dry snow measurements, the cold room was turned off and the door was opened to ambient air, gradually increasing the air 195 temperature in the cold room to room temperature (20 o C). The NIR-HSI images were taken every 1-3 minutes during the warming process and SLF sensor measurements, in each ROI, were taken in between images. For comparison between the two instruments, the mean LWC was calculated over the 10,717 pixels within each ROI, resulting in a 0.4 mm 2 resolution, and was compared against the corresponding SLF sensor measurements at each time step, creating a densely populated comparison dataset spanning a wide range of LWC.

Field
To demonstrate the applicability of the NIR-HSI method for retrieving LWC and re in the field, natural snow was imaged across the vertical wall of a snowpit. Images of the wall were taken at 13:00, at which time there were few clouds, and the air temperature was 210 10 o C. The imager was mounted on to the Resonon Outdoor Field System which includes a tripod mounted rotational stage. The snowpit wall was illuminated with two 500-watt halogen lamps mounted on a tripod, which was placed perpendicular to the wall, similar to the laboratory setup presented in Donahue et al. (2021). For controlled lighting conditions, sun light (direct and diffuse) was blocked by placing an opaque tarp over the top of the snowpit. A detailed schematic of the field setup is shown in Figure 4. For a pixel-by-pixel calibration of the 215 NIR-HSI measurements from radiance to reflectance, a 36% spectrally flat reflectance calibration tarp was hung in front of the snow pit wall, completely covering the field of view of the imager and ROI of the snowpit. The entire vertical profile of the snowpit was not captured in a single image, therefore two images were taken and stitched together capturing only the upper 110 cm of the snowpit at a 6.5 mm 2 pixel resolution. The bottom 40 cm of the snowpit was not captured because perpendicular illumination conditions could not be achieved with the minimum 220 height of the lighting and imager tripods used.
Immediately following imaging, LWC measurements were made with the SLF sensor along a vertical profile at 5 cm increments. For optimal LWC measurements, the SLF sensor requires the dry snow density,  however, the snow was already wet and therefore the dry snow density could not be obtained. Instead, the density of the snow from the adjacent density cut measurement was used, which introduced a small error in the LWC 225 measurement (FPGA Company, 2018) further discussed in Sect. 5.1.

Single scattering
To simulate the optical properties of snow, the single scattering optical properties of constituents (ice, air, water, and impurities), as well as their relative arrangement to one another must be represented. The scattering properties of a 230 single particle are described using three dimensionless optical parameters: (1) the absorption efficiency , (2) scattering efficiency , and (3) asymmetry factor (Bohren and Huffman, 2008). Dry snow particles are often assumed to scatter as a collection of spheres with radius re (e.g., Nolin and Dozier (2000)). With a known re and complex refractive index, Mie scattering theory (Bohren and Huffman, 2008) can be used to calculate the three dimensionless optical parameters. For clean dry snow, this modeling approach is straightforward because only the 235 complex refractive index of ice is needed. For wet snow, on the other hand, the complex refractive index of ice and water is volume-mixed, and there are multiple approaches that can be used to define the arrangement of ice and water relative to each other. The three previously proposed mixing models used in this comparison, defined in Sect.
First, the keff sphere model uses a collection of spheres with re and keff, which was determined using 240 volume-weighted portions of the complex refractive index of ice (kice) and water (kwater).
Second, the coated sphere model used the Coated-Mie scattering code of Mätzler (2002). The radius of the ice core (ri) and thickness of the water coating (tw) was determined by volume-weights of ice ( ) and water ( ).
r e = r + t w Lastly, the interstitial sphere model calculates the single scattering optical properties (i.e., , , and ) of pure 245 ice spheres and water spheres separately, each having the same re. Then, the optical properties were mixed using a volume-weighted average. This method is similar to the keff spheres model; however, the differences are a result of the non-linearity of Mie scattering theory. At shorter wavelengths, the differences in reflectance are small, but at larger wavelengths a notable divergence occurs, shown in Figure 5.

Multiple scattering
To generate a spectral library to match to measured spectra, directional-hemispherical reflectance for each mixing 255 model was simulated using a general-purpose 16-stream plane-parallel discrete ordinates radiative transfer model, DISORT (Stamnes et al., 1988). DISORT allows the user to define optical properties of multiple layers; here, a single optically thick layer was used since the penetration of NIR light into the snowpack is shallow. Optical property inputs for this layer included single scattering albedo, defined as the ratio of the scattering efficiency and extinction efficiency ( + ), and , which were computed from Mie scattering theory. Incoming light can be 260 modeled at multiple user defined zenith angles; here, the zenith angle was set to 0 o to represent nadir lighting.
The output from DISORT was directional-hemispherical reflectance, whereas NIR-HSI measurements are bidirectional reflectance. This is a suitable approach because of the experimental setup: the NIR-HSI measurements were made at nadir illumination and viewing angles and calibrated using a Lambertian white reference target. Under these conditions, snow is nearly Lambertian, allowing for a direction comparison, which would not be the case for 265 non-nadir viewing angles, given that snow heavily favors forward scattering (Dumont et al., 2010).

Retrieving snow properties
To simultaneously retrieve LWC and re, measured spectra in each pixel of the NIR-HSI image was compared to each spectrum in the spectral library to determine the best match by finding the minimum least square residual across 106 bands, ranging from 961 nm to 1472 nm. Hereafter, this method for matching measured and simulated spectra is 270 referred to as the "residual method". An example of (1) measured spectrum from NIR-HSI, (2) retrieved simulated spectra using the interstitial sphere model, and (3) residuals across each band are shown in Figure 6. Once the measured spectrum is best matched to a modeled spectrum, the associated re and LWC are assigned to the pixel, producing separate maps of re and LWC for the imaged area. To reduce impact from sensor noise at the lower limit of the sensor (900 nm), 961 nm was chosen as the starting point for calculating the residuals while still capturing 275 ample spectral data for the left side of the absorption feature centered at 1030 nm ( Figure 6). At longer wavelengths, 1472 nm was chosen as the endpoint for calculating residuals because both ice and water are highly absorptive beyond this point, resulting in a low signal to noise ratio.
Additionally, re was mapped at each timestep using the scaled band area method following Donahue et al. (2021). Briefly, the scaled band area is the area underneath a continuum line spanning an absorption feature and is a 280 shape-based method that is independent of absolute reflectance. Here, the scaled band area was calculated for measured spectra using a predefined start and end point for the continuum line across the absorption feature centered at 1030 nm. The pixel-by-pixel calibration performed here reduced noise at the lower limit of the sensor compared to Donahue et al. (2021), allowing for the defined continuum endpoints to be similar to bands suggested by Nolin and Dozier (2000) (i.e. 961 and 1087 nm). A look up table was populated with "dry snow" (all ice, no water) scaled 285 band areas for simulated re ranging from 30-1500 μm. We note that we used the simulated spectra from the interstitial sphere model, but that the dry snow representations for all mixing models were identical, with variation in spectra introduced only when water was represented. This allowed us to (1) define the starting re for each of the prepared laboratory samples, (2) compare re retrieved from scaled band area to that from the residual method, and (3) assess the suitability of the scaled band area method for grain size retrievals over wet snow.  The full performance comparison across all datasets is summarized in Figure 8, which plots LWC from the 325 SLF sensor against that from each mixing model applied across all samples. To help visualize the difference between experiments, the comparison points are symbolized by different colors as well as marker size that varies with the mean initial dry snow re retrieved using the scaled band area method. Close proximity to the 1:1 line would indicate the best match between NIR-HSI and the SLF sensor. The root mean square error (RMSE) and bias for each sample is summarized in Table 2.

330
For the two samples with the smallest initial re (113 and 130 μm), precipitation particles and decomposed/fragmented particles, it was found that the LWC retrieval method did not perform well using any of the mixing models, with the highest RMSE and bias ( Table 2). The retrieval is near 0% LWC in dry snow conditions, however the LWC retrievals increase rapidly as small amounts of water are introduced. In the coated spheres model, the retrievals reach the LWC limit of the spectral library (25%) when the measured LWC from the SLF sensor was 335 ~10% ( Figure 8B). For keff and interstitial spheres, there is a similar pattern of LWC increasing too fast, although these models do not reach the upper limit of the spectral library. Because none of the models evaluated performed well for small grain sizes, these two samples (113 & 130 μm) were excluded from the calculation of best fit line (red line in Figure 8). This result is discussed further in Sect. 5; however, the exclusion of this data is reasonable because water is not commonly mixed with these types of particles (new snow) in natural environments and therefore not a 340 primary focus for wet snow mapping applications.
For the remaining samples, ranging from 176 -898 μm, the mixing models retrieve LWC values that more closely match the SLF sensor measurements. Visually, both keff and interstitial spheres fall close to the 1:1 line. The keff spheres do have the lowest RMSE and bias for the smallest grains (samples 1-3), but for the remaining medium to large grain size samples, the retrieved LWC consistently has a negative bias and the RMSE is ~2%. Overall, the 345 interstitial spheres have the lowest RMSE (~1%) and bias, which does not trend with grain size, with the best comparison for samples 4-7 and similar values to keff for sample 3. The uncertainty of the model was determined by taking the mean RMSE across samples 3-7, which is 1.4%. This is supported by the dry snow retrieval shown in Figure 7A, where no pixels retrieved greater than 1% LWC. The coated spheres model performed most poorly relative to measurements with a high RMSE and consistently had a positive bias across all samples, although like the 350 other mixing models, the values do fall close to the 1:1 line at low LWC (< 7%). Overall, these experiments and model comparisons show that the interstitial spheres model performed exceptionally well for medium to large grains, and LWC ranges between 0% and 15%, which are the conditions most likely to be found in natural snow covers.

355
Using the residual method, all mixing models retrieved similar grain size values because the grain size retrievals are primarily dependent on the absolute reflectance driven by ice absorption. Here, we present results from the interstitial sphere model because it performed best in the LWC retrieval. For initial, dry snow, conditions the re retrieval using the scaled band area method (re, SBA) had a positive bias relative to the residual method (re, residual), and the bias becomes more positive with increasing grain size (Figure 9). For the remaining, wet snow comparisons, the 360 re, SBA remains relatively constant or increases at low LWC followed by a decrease at high LWC. For Samples 1 and 2, re, SBA remained flat with increasing LWC. For Samples 3-6, re, SBA increased initially with low LWC and then decreases at high LWC. For Sample 7, re, SBA slightly increases before significantly decreasing with increasing LWC.
Whereas the re, residual increases with increasing LWC for all samples, which is expected because the presence of water is known to accelerate snow grain growth (Marsh, 1987). The comparison, further discussed in Sect. 5.2,

365
shows that the scaled band area method is heavily impacted by the presence of water, such that the residual method may be better suited for wet snow.

Field experiment
Although a controlled laboratory environment is ideal for identifying the best suited mixing model, the primary applications for this method would be in situ field studies, motivating the field-based testing of the re and 370 LWC retrieval. The maps of re and LWC generally reflect the profile measurements of LWC and stratigraphy, though at significantly higher detail ( Figure 10). The snowpit was representative of a spring intermountain snowpack undergoing melt, with the density values ranging from ~300 to ~450 3 , and LWC measurements ranging from ~5-15%. The temperature profile, not shown, was isothermal at 0 o C. The general stratigraphy was ice lenses and melt form grains in the upper layers, rounded grains in the central portion, and faceted grains (depth hoar) near 375 the ground. Note that the full pit profile is shown for observations ( Figure 10A and 10D), while the re and LWC retrievals extend across only the upper 110 cm ( Figure 10B and 10C).
The maps were processed using each of the mixing models introduced above, but based on the laboratory findings, only retrievals from the interstitial sphere model are presented here. The SLF sensor measurements were taken along the left side of the ruler (gray stripe or NaNs down the center of snowpit), represented as the dashed red 380 box in Figure 10C. For comparison, LWC was depth averaged in pixels covering the area of the SLF sensor measurements (red line in Figure 10D), in addition to the depth average LWC across the entire width of the snowpit (gray line in Figure 10D). The mean LWC, standard deviation (σ) and number of measurements (n) for the SLF sensor, NIR-HSI along the same profile as the SLF sensor, and the depth average of the entire width of the snowpit are shown in Table 3. Generally, the LWC from the NIR-HSI retrieval had a positive bias compared to the SLF 385 sensor, particularly notable at the top of the snowpack, although the patterns and peaks have similar trends. This discrepancy is attributed to not knowing the dry snow densities needed for the SLF sensor, which would create a negative bias in the measured LWC. The difference between retrievals and measurements is discussed further in section 5.1.
For comparison, when processed with the coated spheres model the LWC reached the model limit at 25% 390 over much of the scene. The keff spheres model did have mean values slightly closer to the SLF sensor, however this was not unexpected given that the laboratory results were biased negative. The magnitude of values, with coated spheres retrieving higher values, and keff spheres retrieving lower values, relative to interstitial spheres mimics the bias present in laboratory results ( Table 2).
The high resolution of the maps shows how stratigraphy influences re and LWC distributions in higher 395 detail than can be captured with standard field-based observations. At the top of the snowpit, it can be discerned that the snow was saturated with water and was pooling on top of the ice lens at 130 cm snow height. The SLF sensor  Water pooling above ice lenses, rather than the ice lens itself having water content, is sensible because there is 400 minimal pore space for water to reside or pass through. Below the ice lens, a few isolated preferential flow paths extend to lower layers, while other features show water concentrated (at ~90 cm, and at ~75 cm) but notably not flowing along the plane of the snowpit wall.

405
The interstitial and keff spheres models performed similarly because the optical properties are volume mixed in both cases, albeit internally versus externally mixed. The external mixing of interstitial water spheres results in a notable shift in the reflectance spectrum at wavelengths ranging from 1300 -1450 nm, when compared to keff spheres ( Figure 5). Since the interstitial sphere model performed the best, this result indicates that the particle size of water in wet snow plays an important role in the simulated reflectance, whereas the particle size of water is not considered 410 in the keff model. Before mixing the optical properties, the particle size of the water sphere was the same size as the ice sphere, which is a reasonable approximation. It would be possible to mix water and ice spheres of differing size, although this approach is computationally expensive, given that the number of possible simulated spectra combinations would approach infinity.
The coated sphere model performed reasonably in the pendular regime, but then overestimated LWC in the 415 funicular regime. The coated sphere model was chosen over the interstitial sphere model by Green et al. (2002), based on visual inspection, considering the bands in the ice absorption feature centered at 1030 cm, which encompasses only part of the distinct shifts between ice and water that are present in the complex refractive index across the NIR (Figure 1). This study used a greater number of NIR bands that span multiple distinguishable shifts between ice and water, which is a more robust approach.

420
For small snow grains of PP (Sample 1) and DF (Sample 2) crystal type, LWC retrievals did not perform well using any of the models. One potential reason being that PP and DF crystal types are complex shapes and reflectance may not be accurately represented using spheres of ice. Using a ray tracing model, Picard et al. (2009) showed that grain shape can influence reflectance. Although it is possible to have wet PP and DF crystal types (e.g., rain on snow), low density dendritic snow crystals are more commonly found at temperatures well below freezing 425 (Judson and Doesken, 2000). It is far more common for wet snow to contain larger rounded grains primarily because the presence of water rapidly increases the rate of snow grain growth, especially through melt and refreeze cycles. Interestingly, for the largest grains, samples 6 and 7, the highest LWC measurements from the SLF sensor are ~10% (Figure 8). These apparent maxima, seen in both measurements and retrievals, is attributed to the large grains having a reduced water holding capacity within the pore space (Yamaguchi et al., 2010). Although the snow sample could contain higher than 10% LWC by volume, water is able to drain below the near surface detection limit 435 of the SLF sensor and NIR-HSI.

Effective grain radius retrieval
The scaled band area method assumes dry snow, but in remote sensing and field applications there is typically no a priori knowledge of snow wetness, thus comparing re, SBA to re, residual allows us to test the validity of this assumption for wet snow. To visualize the difference between the residual and scaled band area methods, an example measured 440 spectrum from a dry and wet snow (12% LWC) sample are shown in Figure 11, along with the corresponding simulated spectrum retrieved using the two methods. For the dry snow spectrum, re, SBA is larger than re, residual, which is a consistent trend across all samples (Figure 9 best fit line). This is attributed to the scaled band area method using a fixed wavelength range spanning the area of the absorption feature centered at 1030 nm (961 to 1087 nm), which makes it sensitive to small changes in the shape of the measured absorption feature, and location of the continuum 445 line shoulders (shown as gray lines in Figure 11B). Whereas the residual method finds the best fit spectrum over all NIR wavelengths ranging from 961 nm to 1472 nm. In the example shown ( Figure 11B), the left shoulder of the matched simulated spectrum using the residual method is slightly below the measured spectrum, resulting in the scaled band area of the measured spectrum (2.92) being higher than that of the matched simulated spectra (2.46).
In the wet snow case, the absorption features centered at 1030 nm shifts to shorter wavelengths and 450 broadens. Similar to the comparison for dry snow, the fixed wavelength range and continuum line fails to fully capture the wet snow absorption feature, notably resulting in a reduction of the scaled band area. This result is shown in the spectral reflectance example ( Figure 11C) and is responsible for the decreasing re, SBA at high LWC  Figure 11C), resulting in a smaller scaled band area. This example also shows the broadening of the feature, relative to the shape of ice absorption alone. The shape of the wet snow feature is better represented by the matched simulated spectrum using the residual method. It is possible that this has implications for previous studies that have applied the scaled band area method over potentially wet snow (e.g., Skiles and Painter (2017)).

460
All the mixing models examined in this study are approximate representations of the relative arrangement of ice and water in wet snow. The spherical particle approximation used in this study to represent wet snow is a reasonable approach because ice grains in the presence of water tend to be rounded. The arrangement of ice and water, on the other hand, is dependent on the level of water saturation, therefore using a single mixing model to determine the LWC across a large range of water saturations results in inherent uncertainty. Since no a priori knowledge of snow 465 wetness is known when taking NIR-HSI measurements, the goal of this research only aims to find the modeling approach that has the best retrieval of LWC across ranges commonly found in natural environments, when compared to an established dielectric method-based instrument.
Dielectric instruments, including the SLF sensor used here, have their own uncertainties. These types of instruments rely on an independent snow density measurement to calculate LWC. In the case of the SLF sensor, the 470 "dry snow" density, which describes the fraction of ice in a wet snow volume, is required to calculate a LWC value (FPGA Company, 2018). In the laboratory experiments, this was not an issue because all snow samples started dry, and the dry snow density was used for all subsequent LWC measurements during melt. Conversely, in the field experiment, the snow was already wet and there was no way of isolating a dry snow density. One solution offered by FPGA Company (2018) is to measure density in the morning if the snow has refrozen, but this was not possible.

475
Although not quantifiable here, the error introduced is assumed to be small since the measured LWC is dominated by water content rather than snow density (FPGA Company, 2018). Nevertheless, this error could attribute to the discrepancy in LWC retrieved from NIR-HSI and the SLF sensor seen in Figure 10D. A dry snow density, if known, would be lower than the measured wet snow density, resulting in a slightly higher LWC measurement and closer match to the NIR-HSI retrievals. Although there is inherent uncertainty associated with the SLF sensor, it was found 480 to be the most suitable LWC measurement instrument for comparison to LWC retrieved from NIR reflectance because it is non-destructive to the snow surface, measures a flat surface, and has minimal penetration into snow, similar to NIR light.
Lastly, there is also uncertainty in the measurement of absolute spectral reflectance, the accuracy of which is important for using the residual method. To minimize this uncertainty in the laboratory, spectral measurements were 485 taken with consistent lighting conditions and at close proximity, resulting in near perfect conditions, which was ideal for the comparison study. Additionally, the residual method benefits from the high signal-to-noise ratio (1885) and spectral resolution (4.9 nm) of the instrument used here. For application of this method at the airborne or satellite scales, spectral measurements contain more noise than those in the laboratory and require atmospheric and https://doi.org/10.5194/tc-2021-247 Preprint. Discussion started: 12 August 2021 c Author(s) 2021. CC BY 4.0 License.
topographic correction, introducing additional uncertainty into absolute reflectance values. Although not within the 490 scope of this study, instruments at these scales also need to account for water vapor in the air which is discussed further in Green et al. (2006). Finally, due to the relatively minor shift in NIR reflectance, this approach is probably not suitable for mixed pixels (not pure snow), which become more common as spatial coverage and pixel sizes increase.

495
The results in this study show the externally mixed interstitial spheres model performs best when compared to a dielectric LWC measurement instrument, relative to the previously proposed keff spheres (Hyvarinen and Lammasniemi, 1987) and coated spheres model (Green et al., 2002). It was also found that for the smallest grains (i.e., new and decomposing precipitation particles) none of the models investigated compared well to the SLF sensor. For low LWC (<~7%), all of retrievals compare well to measurements, but at higher LWC the keff and coated 500 spheres were biased positive and negative, respectively. Overall, and across the widest range of initial grain types, re (162 -859 μm) and LWC (0-17.2%), the interstitial sphere model performed the best, with ~1% uncertainty.
For the re comparison, between the two optically based residual and scaled band area methods, it was found that the scaled band area retrieval had a positive bias compared to the residual method. This bias was lowest for small grains and increased with grain size. For wet snow, the scaled band area method was impacted by the presence of 505 water due to the shift in the absorption feature to shorter wavelengths, resulting in decreasing re at high LWC.
Because scaled band area implicitly assumes dry snow, it is based on ice absorption alone, caution is encouraged when applying this method without knowing that the snow is dry.
The field application of the NIR-HSI method produced maps that reflect the general understanding of what a snapshot of a snowpit would look like during snowmelt progression, but mapped the stratigraphy of snow properties 510 (i.e., LWC and grain size) at much higher detail relative to standard profile-based observations. The retrieved LWC was found to be slightly higher than that measured by the SLF sensor, which was attributed to both the inability to determine the dry snow density and the high level of detail in the maps that could not be captured by the volume averaged dielectric sensor.
For future work and deployment of NIR-HSI, viewing re and LWC spatial distributions in such unprecedented 515 detail has implications for better understanding water movement patterns through snow, how re and LWC relate to other snow properties and topographic characteristics, and for interpreting microwave remote sensing retrievals over wet snow. There is also broader relevance for the assessment and development of snow property retrievals from measured spectral reflectance with upcoming satellite imaging spectrometer missions. These include the Surface Biology and Geology (SBG) imaging spectrometer mission and the Copernicus Hyperspectral Imaging Mission 520 (CHIME). Although algorithm suites have been developed to retrieve snow properties from airborne imaging spectroscopy (Painter et al., 2013), LWC is not a standard part of the retrieval, and has only been demonstrated in case studies (Green et al., 2002).

Data availability
Data will be available upon request from lead author during review and will be placed into a data repository upon 525 acceptance.

Author contribution
CD conceptualized the study, collected laboratory and field data, conducted DISORT model simulations, and analyzed the results. MS provided the hyperspectral imager, assisted with data collection, and advised CD during conceptualization and DISORT model simulations. KH acquired funding for this research, was responsible for 530 project administration, and supervised CD throughout the study. CD wrote the manuscript with contributions from all coauthors.

Competing interests
The authors declare that they have no conflicts of interest.