Time-resolved measurements of the densities of individual frozen hydrometeors and fresh snowfall

It is a challenge to obtain accurate measurements of the microphysical properties of delicate, structurally complex, frozen, and semi-frozen hydrometeors. We present a new technique for the real-time measurement of the density of freshly fallen individual snowflakes. A new thermal-imaging instrument, the Differential Emissivity Imaging Disdrometer (DEID), has been shown through laboratory and field experiments to be capable of providing accurate estimates of individual snowflake and bulk snow hydrometeor density (which can be interpreted as the snow-to-liquid ratio or SLR). The method exploits the rate of heat transfer during the melting of a hydrometeor on a heated metal plate, which is a function of the temperature difference between the hotplate surface and the top of the hydrometeor. The product of the melting speed and melting time yields an effective particle thickness normal to the hotplate surface, which can then be used in combination with the particle mass and area on the plate to determine a particle density. Uncertainties in estimates of particle density are approximately 4 % based on calibrations with laboratory-produced particles made from water and frozen solutions of salt and water and field comparisons with both high-resolution imagery of falling snow and traditional snowpack density measurements obtained at 12 h intervals. For 17 storms, individual particle densities vary from 19 to 495 kg m⁻³, and storm mean snow densities vary from 40 to 100 kg m⁻³. We observe probability distribution functions for hydrometeor density that are nearly Gaussian with kurtosis of ≈ 3 and skewness of ≈ 0.01.

1 Introduction

Frozen and semi-frozen hydrometeors have a very wide range of porosities (Dunnavan et al., 2019). Determining their particulate densities and bulk snow-to-liquid ratios (SLRs) once fallen on the ground is important to a wide range of fields, including hydrology (Rango and Martinec, 1995; Sturm et al., 2010), climatology (Dickinson, 1983), remote sensing at wavelengths ranging from the visible to the microwave (Kendra et al., 1994; Kokhanovsky and Zege, 2004; Gergely et al., 2010), and the parameterization of snowflake fall speeds in weather and climate models (Rutledge and Hobbs, 1984; Hong et al., 2004; Fovell and Su, 2007; Alcott and Steenburgh, 2010; Finlon et al., 2019). Because hydrometeor porosity is invisible to most imaging techniques, obtaining accurate snowflake density estimates has proven to be a significant challenge where even the best estimates have required the use of sophisticated field programs using multiple instruments (Tiira et al., 2016).

Improvements to weather prediction are currently hampered by an inability to assimilate information about ongoing variability in frozen and semi-frozen precipitation particles (Rasmussen et al., 2011). Avalanche forecasting in mountainous regions depends, in part, on knowledge of the vertical density structure of freshly fallen snow (Morrison et al., 2023), a parameter that is typically measured at sparse intervals (Schweizer et al., 2011; Proksch et al., 2016) using techniques such as micro-computed tomography (μCT) or, more typically, manual gravimetric methods (Proksch et al., 2016).

In our previous work, we showed that a new thermal-imaging instrument, the Differential Emissivity Imaging Disdrometer (DEID), can be used to measure individual hydrometeor density based on the first automated direct measurements of particle mass in combination with estimates of the spherical-particle-equivalent effective diameter or by using concurrent photographic imagery of the morphological characteristics of hydrometeors as they fall (Singh et al., 2021; Rees et al., 2021). While the spherical-particle approach offers the advantage of simplicity, it was found to lead to snowflake density estimates that were significantly biased low relative to a method that required an added camera system, likely because snowflakes are not in fact spheres. Here, we describe a new method for estimating particle-by-particle frozen-hydrometeor density that, like the spherical-particle method, uses only the DEID to measure mass but instead infers particle volume from DEID measurements of melting time and particle area and estimates of the rate of heat transfer from the hotplate to the hydrometeor to obtain a melting speed (MS).

2 DEID measurement techniques for obtaining hydrometeor mass and density

The DEID consists of an infrared camera pointed at the surface of a low-emissivity aluminum hotplate. To quantify hydrometeor area on the hotplate, the DEID makes use of the contrasting thermal emissivities of water (ε>0.95) and aluminum (ε<0.1) at the same temperature. Owing to the high difference in emissivity, melted hydrometeors with nearly the same thermodynamic temperature as the heated plate have strongly contrasting radiative temperatures such that droplets on the heated plate can be easily discriminated using a thermal camera. The hotplate surface is roughened, which prevents displacement of melted snowflakes at high wind speeds, as demonstrated in wind-tunnel experiments with wind speeds varying from 2 to 12 m s⁻¹.

2.1 Particle mass measurement

The DEID methodology for obtaining the mass of a hydrometeor particle has previously been described by Singh et al. (2021), Rees et al. (2021), Rees and Garrett (2021), and Morrison et al. (2023). Here, we present a concise summary including recent modifications to the measurement methodology. Briefly, the mass of individual hydrometeors is obtained by considering the contact area of each hydrometeor on a heated metal plate and the temperature difference between the plate and the surface of the melted liquid particle, which is integrated over time from the point of first impact of the particle onto the plate surface up to the point of its complete evaporation.

Figure 1 For the methodology presented here, a control volume is defined to wrap around the hydrometeor. (a) Side view of control volume. Schematic of the heat-transfer process from the DEID's aluminum hotplate to a solid hydrometer along the z axis during melting. R is the radius of the hemispherical ice particle on the hotplate, h is the effective thickness of the hemisphere that is $\mathrm{2}R/\mathrm{3}$, and m is the mass of the ice particle. (b) Top view of the control volume. A_c is the contact area of the ice particle in the x–y plane. (c) Schematic illustrating the method for calculating the contact area of a frozen hydrometeor with an arbitrary geometry on the hotplate. The individual rectangles represent pixels as viewed with the thermal camera. (d) Schematic illustrating the method for calculating the volume of a hydrometeor using the contact area and height of a frozen hydrometeor with an arbitrary geometry on the hotplate. $h=(\mathrm{1}/N){\sum }_i{\sum }_j{h}_{ij}$, where h_ij is the height of a frozen hydrometeor normal to the hotplate that is associated with the ijth pixel at t=0, and N is the total number of pixels associated with a frozen hydrometer at t=0.

Specifically, the mass m of an individual hydrometeor is obtained by applying conservation of energy to a control volume surrounding each hydrometeor on the hotplate (see Fig. 1). The heat gained by a snowflake from the heated plate is assumed to be equivalent to the heat required to increase the snowflake's internal energy and the heat lost during melting and evaporation, which may be described as follows:

$$\begin{array}{}\text{(1)}& Q_{\mathrm{in}}-Q_{\mathrm{out}}=\mathrm{\Delta }{Q}_{\mathrm{st}}.\end{array}$$

Here, Q_in is the conductive heat gained by a hydrometer from the hotplate, Q_out is the heat loss from a hydrometeor to the surroundings, and ΔQ_st is the energy stored in a hydrometer during melting and evaporation. Neglecting convection and radiation losses from the hydrometer (which was shown to be a good approximation in Singh et al., 2021), Eq. (1) can be written as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \underset{t_{\mathrm{pc}}}{\int }\underset{A_{\mathrm{c}}}{\int }k(\partial T/\partial z)\mathrm{d}{A}_{\mathrm{c}}\mathrm{d}t-\underset{\mathcal{V}}{\int }\mathit{\rho }{L}_{\mathrm{eqv}}\mathrm{d}\mathcal{V}\\ & =\underset{t_{\mathrm{pc}}}{\int }\underset{\mathcal{V}}{\int }\mathit{\rho }{C}_p(\partial T/\partial t)d\mathcal{V}\mathrm{d}t.\end{array}\end{array}$$

Here, k is the thermal conductivity of the aluminum plate, A_c is the contact area of the hydrometeor on the hotplate, t_pc is the time it takes to melt or evaporate a droplet (phase change), $\partial T/\partial z$ is the temperature gradient related to conduction from the hotplate into the hydrometeor, ρ is the density of the hydrometeor, 𝒱 is the volume of the hydrometeor, $L_{\mathrm{eqv}}=L_{\mathrm{f}}+L_{\mathrm{v}}$ is the combined latent heat of fusion (L_f) and vaporization (L_v), C_p is the specific heat of the hydrometeor, and z is the vertical direction normal to the hotplate. Using $m={\int }_{\mathcal{V}}\mathit{\rho }d\mathcal{V}$, $\partial T/\partial z\approx \mathrm{\Delta }T(x,y,t)/\mathrm{\Delta }z$ and by separating the energy stored in the hydrometeor during melting and evaporation, Eq. (2) can be modified as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \underset{t_{\mathrm{pc}}}{\int }\underset{A_{\mathrm{c}}}{\int }(k/\mathrm{\Delta }z{)}_{\mathrm{eff}}\mathrm{\Delta }T(x,y,t)\mathrm{d}{A}_{\mathrm{c}}\mathrm{d}t-m(L_{\mathrm{f}}+L_{\mathrm{v}})\\ & =m{C}_{\mathrm{ice}}\left(T_{\mathrm{0}}-T_{\mathrm{ice}}\right)+m{C}_{\mathrm{liq}}\left(T_{\mathrm{p}}-T_{\mathrm{0}}\right).\end{array}\end{array}$$

Here, C_ice is the specific heat of ice; T_p is the hotplate surface temperature during melting and evaporation (it is constant with time); T_ice is the initial temperature of a frozen hydrometer; T₀=0 °C; $L_{\mathrm{f}}=\mathrm{3.34}\times {\mathrm{10}}^{\mathrm{5}}$ J kg⁻¹ is the latent heat of fusion of water; $L_{\mathrm{v}}=\mathrm{2.32}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹ is the latent heat of vaporization of water (see Appendix A1); $C_{\mathrm{liq}}=\mathrm{4.18}\times {\mathrm{10}}^{\mathrm{3}}$ J kg${}^{-\mathrm{1}}{K}^{-\mathrm{1}}$ is the specific heat of liquid water; and $(k/\mathrm{\Delta }z{)}_{\mathrm{eff}}=\mathit{\kappa }$ is an empirical device-specific calibration coefficient related to the amount of heat that passes through the metal plate into individual hydrometeors per unit of time through a unit area with a temperature gradient of 1 °, determined to be $\mathrm{7.01}\pm \mathrm{0.01}\times {\mathrm{10}}^{\mathrm{3}}$ W m⁻² K⁻¹ (see Singh et al., 2021, for details), which is independent of particle size and environmental conditions. In practice, to numerically calculate the mass of a hydrometeor using images from a thermal camera on a “pixel-by-pixel” basis, we suppose that n_x is the total number of pixels making up the hydrometeor in the x direction at time t, and n_y is the total number of pixels associated with the hydrometeor in the y direction at time t. The total number of pixels at time t is N(t)=n_x(t)n_y(t). Equation (3) can be modified and rearranged for a pixel-by-pixel implementation of Eq. (3) as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& m{C}_{\mathrm{ice}}(T_{\mathrm{0}}-T_{\mathrm{ice}})+m{L}_{\mathrm{f}}+m{C}_{\mathrm{liq}}\left(T_{\mathrm{p}}-T_{\mathrm{0}}\right)+m{L}_{\mathrm{v}}\\ & =\mathit{\kappa }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{\mathrm{evap}}}{\int }}\sum _i^{n_x}\sum _j^{n_y}\left(T_{\mathrm{p}}-T_{ij}\left(t\right)\right)A_{ij}\left(t\right)\mathrm{d}t.\end{array}\end{array}$$

Here, Δt_evap is the time required to melt and evaporate a hydrometeor, A_ij(t) is the area of the ijth pixel at time t, and T_ij(t) is the temperature of the ijth pixel at time t. Assuming all pixels have the same area (e.g., A₁), then ${\sum }_i^{n_x}{\sum }_j^{n_y}{A}_{ij}\left(t\right)=n_x{n}_y{A}_{\mathrm{1}}\left(t\right)$ is the total contact area of the hydrometeor at time t (i.e., n_xn_yA₁(t)=A_c(t)). Substituting A_c(t) into Eq. (4) yields

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& m{C}_{\mathrm{ice}}\left(T_{\mathrm{0}}-T_{\mathrm{ice}}\right)+m{L}_{\mathrm{f}}+m{C}_{\mathrm{liq}}\left(T_{\mathrm{p}}-T_{\mathrm{0}}\right)+m{L}_{\mathrm{v}}\\ & =\mathit{\kappa }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{\mathrm{evap}}}{\int }}{A}_{\mathrm{c}}\left(t\right)(T_{\mathrm{p}}-{\displaystyle \frac{\mathrm{1}}{n_x{n}_y}}\sum _i^{n_x}\sum _i^{n_x}{T}_{ij}(t\left)\right)\mathrm{d}t.\end{array}\end{array}$$

Here, $\frac{\mathrm{1}}{n_x{n}_y}{\sum }_i^{n_x}{\sum }_i^{n_x}{T}_{ij}\left(t\right)=T_{\mathrm{h}}\left(t\right)$. T_h(t) is the spatial mean temperature over all frozen hydrometeors and/or water droplet pixels at time t on the hotplate. Substituting the spatial mean temperature T_h(t) into Eq. (5) yields

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& m{C}_{\mathrm{ice}}\left(T_{\mathrm{0}}-T_{\mathrm{ice}}\right)+m{L}_{\mathrm{f}}+m{C}_{\mathrm{liq}}\left(T_{\mathrm{p}}-T_{\mathrm{0}}\right)+m{L}_{\mathrm{v}}\\ & =\mathit{\kappa }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{\mathrm{evap}}}{\int }}{A}_{\mathrm{c}}\left(t\right)\left(T_{\mathrm{p}}-T_{\mathrm{h}}\left(t\right)\right)\mathrm{d}t.\end{array}\end{array}$$

In Eq. (6), at any time t during melting and evaporation, the temperature of some portions of the hydrometeor area is less than 0 and some temperatures are greater than 0. The area of a hydrometer at any time t can be written as $A_{\mathrm{c}}\left(t\right)=A_{\mathrm{ice}}\left(t\right)+A_{\mathrm{liq}}\left(t\right)$, where A_ice(t) is the contact area of the hydrometeor fraction on the hotplate with a temperature less than or equal to 0 at time t (the sum of all pixels with temperatures less than or equal to 0), and A_liq(t) is the contact area of hydrometeors on the hotplate with a temperature greater than 0 at time t (sum of all pixels with a temperature greater than 0). After substituting A(t) into Eq. (6), Eq. (6) may be re-written as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& m{C}_{\mathrm{ice}}\left(T_{\mathrm{0}}-T_{\mathrm{ice}}\right)+m{L}_{\mathrm{f}}+m{C}_{\mathrm{liq}}\left(T_{\mathrm{p}}-T_{\mathrm{0}}\right)+m{L}_{\mathrm{v}}\\ & =\mathit{\kappa }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{\mathrm{melt}}}{\int }}\left(T_{\mathrm{p}}-T_{\mathrm{ice}}\left(t\right)\right)A_{\mathrm{ice}}\left(t\right)\mathrm{d}t\\ & +\mathit{\kappa }\underset{t_{\mathrm{0}}}{\overset{\mathrm{\Delta }{t}_{\mathrm{evap}}}{\int }}\left(T_{\mathrm{p}}-T_{\mathrm{liq}}\left(t\right)\right)A_{\mathrm{liq}}\left(t\right)\mathrm{d}t.\end{array}\end{array}$$

Here, t=0 corresponds to when a frozen hydrometeor hits the hotplate, and t=t₀ is when the thermal camera sees the liquid portion for the first time. t₀ is the time lag between the melting and evaporation start times, which is about 0.1 s for the laboratory ice particles tested and negligible for snowflakes observed using our typical thermal-camera recording frame rates. T_ice(t) is the spatial mean temperature of all pixels with a temperature less than or equal to 0, and T_liq(t) is the spatial mean temperature of all pixels with a temperature greater than 0. Note that in post-processing data, the thermal camera can be adjusted to selectively “see” particles on the hotplate in specific temperature ranges (see Sect. 3 and Fig. E1). To more easily evaluate Eq. (7), the camera is set to only see hydrometeors on the hotplate after melting (details are given in Fig. E1). We assume that heat transfer through A_ice exclusively leads to increases in ice temperature and melting, and heat transfer through A_liq does not go into these ice portions. We justify this since the temperature gradients between the ice and water are much smaller than the temperature gradients between the plate and hydrometeor. Furthermore, the thermal conductivity of aluminum is much higher than that of ice or water. Hence, A_ice(t)=0 for T_h(t)<0 °C and $\mathit{\kappa }{\int }_{\mathrm{0}}^{\mathrm{\Delta }{t}_{\mathrm{melt}}}(T_{\mathrm{p}}-T_{\mathrm{h}}(t\left)\right)A_{\mathrm{ice}}\left(t\right)\mathrm{d}t=\mathrm{0}$, which is equal to $m{C}_{\mathrm{ice}}(T_{\mathrm{0}}-T_{\mathrm{ice}})+m{L}_{\mathrm{f}}$. This trick allows us to remove both terms from Eq. (7), which yields

$$\begin{array}{}\text{(8)}& m={\displaystyle \frac{\mathit{\kappa }}{C_{\mathrm{liq}}{T}_{\mathrm{p}}+L_{\mathrm{v}}}}\underset{t_{\mathrm{0}}}{\overset{\mathrm{\Delta }{t}_{\mathrm{evap}}}{\int }}\left(T_{\mathrm{p}}-T_{\mathrm{liq}}\left(t\right)\right)A_{\mathrm{liq}}\left(t\right)\mathrm{d}t.\end{array}$$

Here, the lower integration bound is reset so that t₀=0 such that Eq. (8) is evaluated for the liquid phase only, and the time integral is evaluated using the trapezoidal rule. T_liq(t) is the spatial mean water droplet temperature at time t. Note that the initial and final temperatures of all liquid droplets during evaporation are T₀=0 °C and T_p, respectively, and more than ≈ 98 % of the mass evaporates at the highest temperature (≈T_p). This is illustrated in the time series of the temperature of an evaporating hydrometeor given in Fig. A1. Note that we use the subscripts ice and s to denote frozen hydrometeors (i.e., ice and snowflake, respectively).

Mass estimates were shown in wind-tunnel calibrations to be nearly independent of environmental conditions, including wind speed, relative humidity, and ambient temperature (Singh et al., 2021). Specifically, wind-tunnel experiments with the DEID showed less than 4 % variability in mass measurements of hydrometeors for a wide range of wind speeds, relative humidities, and air temperatures (Singh et al., 2021). The reason for the low sensitivity to environmental conditions is that the DEID directly measures the energy required to melt and evaporate a droplet, mL_eqv. For example, the heat-transfer rate to a droplet is dependent on parameters such as wind speed through the Reynolds number (Kosky et al., 2013) and the temperature. However, while higher winds may accelerate heat transfer, they also diminish the time for completing evaporation. Because the product of the heat-transfer rate and evaporation time determines particle mass, winds play a minor role in the calculation of mass.

2.2 Particle density

Obtaining frozen-hydrometeor density from hydrometeor mass requires an estimate of the particle volume while it is in its frozen state. While the DEID can provide an accurate estimate of snowflake mass m and initial snowflake contact area after it impacts the hotplate A_c, it cannot provide a direct measure of a particle's effective thickness in the direction normal to the hotplate h as illustrated in Fig. 1. In its place, we have developed a method for estimating h based on a “melting speed” v_melt such that h=v_meltΔt_melt, where Δt_melt is the time required to melt an individual snowflake. Using these substitutions, the density of a frozen hydrometeor can be written as

$$\begin{array}{}\text{(9)}& {\mathit{\rho }}_{\mathrm{MS}}={\displaystyle \frac{m}{{\mathcal{V}}_{\mathrm{s}}}}={\displaystyle \frac{m}{A_{\mathrm{c}}{v}_{\mathrm{melt}}\mathrm{\Delta }{t}_{\mathrm{melt}}}},\end{array}$$

where ρ_MS indicates the density computed using the melting speed method, and 𝒱_s is the volume of a snowflake estimated as A_cv_meltΔt_melt.

We propose a method to measure v_melt as a function of the temperature difference across a hydrometeor (ΔT_melt) – and hence the heat-transfer rate – as illustrated in Fig. 1. During the time it takes to melt a snowflake, Δt_melt, a hydrometeor receives a quantity of energy equal to $m{L}_{\mathrm{ff}}=m{C}_{\mathrm{ice}}(T_{\mathrm{0}}-T_{\mathrm{ice}})+m{L}_{\mathrm{f}}$, i.e., the sum of the internal energy of a frozen hydrometeor and its latent heat of fusion (energy received by the ice to increase temperature and melt completely) from the hotplate, which is independent of the density of a frozen hydrometeor. v_melt is associated with the conductive heat flux (κΔT_melt) from the hotplate to the frozen hydrometeors during melting. We hypothesize that v_melt is a function of the temperature difference across a hydrometeor (ΔT_melt), and we may write v_melt as

$$\begin{array}{}\text{(10)}& v_{\mathrm{melt}}=c_{\mathrm{melt}}\mathrm{\Delta }{T}_{\mathrm{melt}},\end{array}$$

where $\mathrm{\Delta }{T}_{\mathrm{melt}}=\stackrel{\mathrm{‾}}{T_{\mathrm{p}}-T_{\mathrm{s}}\left(t\right)}$ (the overbar represents a temporal mean) during the melting process, T_s(t) is the spatial mean of the surface temperature of the frozen portion of the particle during melting, and c_melt is a constant determined experimentally (see Sect. 4.1).

Now, if v_melt from Eq. (10) is substituted into Eq. (9), the MS density equation can now be written as

$$\begin{array}{}\text{(11)}& {\mathit{\rho }}_{\mathrm{MS}}={\displaystyle \frac{m}{c_{\mathrm{melt}}{A}_{\mathrm{c}}\mathrm{\Delta }{T}_{\mathrm{melt}}\mathrm{\Delta }{t}_{\mathrm{melt}}}}.\end{array}$$

The melting parameter Δt_melt is quite short for low-density snowflakes and hence requires high-frequency recording of thermal images, resulting in the generation of a tremendous amount of data, which is not convenient for field experiments. In field experiments, the DEID measures Δt_evap, which is much longer than Δt_melt. Fortunately, a relation between Δt_melt and Δt_evap can be derived easily (see Appendix B). By estimating the average conductive heat-transfer rate during the melting and evaporation processes, we may substitute $\mathrm{\Delta }{T}_{\mathrm{melt}}\mathrm{\Delta }{t}_{\mathrm{melt}}\approx (L_{\mathrm{ff}}/L_{\mathrm{vv}})\left(\mathrm{\Delta }{T}_{\mathrm{evap}}\mathrm{\Delta }{t}_{\mathrm{evap}}\right)$ into Eq. (11), which yields

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_{\mathrm{MS}}=c{\displaystyle \frac{m}{A_{\mathrm{c}}\mathrm{\Delta }{T}_{\mathrm{evap}}\mathrm{\Delta }{t}_{\mathrm{evap}}}}.\end{array}$$

Here, $\mathrm{\Delta }{T}_{\mathrm{evap}}=\stackrel{\mathrm{‾}}{T_{\mathrm{p}}-T_{\mathrm{liq}}\left(t\right)}$ during the evaporation process, and the constant, c, is given by

$$\begin{array}{}\text{(13)}& c=\left(L_{\mathrm{vv}}\right)/\left(L_{\mathrm{ff}}{c}_{\mathrm{melt}}\right)\end{array}$$

in units of kelvin per second per meter (K s m⁻¹). The constant, c, is derived from a combination of thermodynamic and laboratory calibration constants that must be determined experimentally (see Sect. 4.1).

In practice, Eq. (11) is evaluated following the methodology shown in Fig. 1b. Snowflakes with complex shapes do not have simple height relationships like $h=\mathrm{2}R/\mathrm{3}$, as shown in Fig. 1a. Hence, a method is required to determine the height and volume of each pixel within a snowflake. If A_ij is taken as the area of the ijth pixel at t=0 and h_ij is assumed to be the height of a frozen hydrometeor normal to the hotplate that is associated with the ijth pixel at t=0, the effective thickness of frozen hydrometer h can be written as $h=(\mathrm{1}/N){\sum }_i{\sum }_j{h}_{ij}$. Here, N is the total number of pixels associated with a frozen hydrometer; h_ij can be estimated using v_meltΔt_ij, where $v_{\mathrm{melt}}=c_{\mathrm{melt}\mathrm{\Delta }{T}_{ij}}$; and 𝒱_ij is the volume of the ijth pixel. We may then write Eq. (14) as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathcal{V}}_{\mathrm{s}}& =\sum _i\sum _j{\mathcal{V}}_{ij}=\sum _i\sum _j{A}_{ij}{h}_{ij}\\ & =c_{\mathrm{melt}}\sum _i\sum _j{A}_{ij}\mathrm{\Delta }{T}_{ij}\mathrm{\Delta }{t}_{ij},\end{array}\end{array}$$

which can subsequently be estimated using the MS method on a pixel-by-pixel basis. Here, ΔT_ij is the temporal mean temperature difference across the ijth pixel, Δt_ij is the melting time of the ijth pixel, and c_melt is the calibration constant of the melting velocity. This study used Eq. (14) to calibrate laboratory ice particles and compare snowflake habits. For field observations, we make the following simplification to determine the volume:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathcal{V}}_{\mathrm{s}}& =c_{\mathrm{melt}}{A}_{\mathrm{c}}{\displaystyle \frac{\mathrm{1}}{N}}\sum _i\sum _j\mathrm{\Delta }{T}_{ij}\mathrm{\Delta }{t}_{ij}\\ & \approx c_{\mathrm{melt}}{A}_{\mathrm{c}}\mathrm{\Delta }{t}_{\mathrm{melt}}{\displaystyle \frac{\mathrm{1}}{N}}\sum _i\sum _j\mathrm{\Delta }{T}_{ij}\\ & =c_{\mathrm{melt}}{A}_{\mathrm{c}}\mathrm{\Delta }{T}_{\mathrm{melt}}\mathrm{\Delta }{t}_{\mathrm{melt}}={\displaystyle \frac{\mathrm{1}}{c}}{A}_{\mathrm{c}}\mathrm{\Delta }{T}_{\mathrm{evap}}\mathrm{\Delta }{t}_{\mathrm{evap}}.\end{array}\end{array}$$

Here, A_c=NA_ij (assuming all pixels have the same area), and $\mathrm{\Delta }{T}_{\mathrm{melt}}=\frac{\mathrm{1}}{N}{\sum }_i{\sum }_j\mathrm{\Delta }{T}_{ij}$ is the spatial and temporal mean temperature difference across a frozen hydrometer during melting. Based on experimental tests, we assume Δt_melt≈Δt_ij, which is the melting time of a frozen hydrometer, and $c=\left(L_{\mathrm{vv}}\right)/\left(L_{\mathrm{ff}}{c}_{\mathrm{melt}}\right)$​​​​​​​ (see Appendix B).

Note that Δt_melt is impacted by variability in environmental conditions. A sample time series of temperature and hydrometeor area during melting and evaporation is shown in Fig. 2. During the melting process, the area of a particle that is in its frozen state decreases to 0 from a maximum immediately after having fallen on the plate. At the same time, the observed liquid component of the hydrometeor increases to a maximum before abruptly disappearing. The sum of these two areas is nearly constant, at least accounting for inevitable uncertainties in the binary thresholding associated with discriminating the hydrometeor from its background. Notably, the sum is also equal to the initial area of the frozen particle prior to its melting. Figure 2 shows how the area of the ice particle prior to melting is similar to the maximum area of the liquid droplet, showing that the area of solid hydrometeors is preserved after melting.

Figure 2 Observed melting and evaporation of a laboratory-made ice particle on the DEID hotplate. (a) Time series of the area of the ice particle, A_ice(t) (dashed black line), and the liquid, A_liq(t) (solid black line). (b) The minimum surface temperature of the ice particle (dashed black line) and the temperature of the liquid, T_liq(t) (solid black line), immediately after being placed on the hotplate. The horizontal dotted line represents a temperature partition in the melting and evaporation process. The shaded region denotes the period of melting. Details on the manufacturer of the laboratory-made ice particles and experiments are provided in Sect. 3.1. The lag between the melting and evaporation start times is 0.1 s, which is equal to the thermal-camera sampling period (sampling rate of 10 Hz).

2.3 Use of the DEID to determine bulk-snowpack-derived quantities

In additional to individual hydrometeor mass and density measurements, the DEID can be used to provide useful bulk snowpack quantities. Precipitation intensity or snow water equivalent rate of precipitation, $\dot{\mathrm{SWE}}$ (in mm h⁻¹), can be estimated from the cumulative particle mass measured by the DEID over a given time period (Δt_res) as

$$\begin{array}{}\text{(16)}& \dot{\mathrm{SWE}}=k{\displaystyle \frac{\mathrm{\Delta }m}{{\mathit{\rho }}_{\mathrm{w}}{A}_{\mathrm{hp}}\mathrm{\Delta }{t}_{\mathrm{res}}}},\end{array}$$

where k is a conversion factor from meters per second to millimeters per hour (3.6 × 10⁶ mm h⁻¹ m⁻¹ s); Δt_res is the sampling time (h); Δm (kg) is the total hydrometeor mass that falls on the hotplate at a given time, where the individual mass of the hydrometeor is estimated using Eq. (8); ρ_w (kg m⁻³) is the bulk density of water; and A_hp (m²) is a rectangular sampling area on the hotplate that captures all hydrometeors. The accumulated SWE (mm) can be calculated over a given time interval Δt_res (h) as $\mathrm{SWE}=\dot{\mathrm{SWE}}\times \mathrm{\Delta }{t}_{\mathrm{res}}$.

The average density (${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}$) of the snowflakes can be estimated using the DEID from the ratio of the cumulative measured mass and the total volume of all snowflakes sampled in a given time interval (Δt_res),

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{m}_i}{{\sum }_{i=\mathrm{1}}^N{m}_i/{\mathit{\rho }}_{\mathrm{MS},i}}},\end{array}$$

where m_i (kg) is the mass of the ith snowflake, ρ_MS,i (kg m⁻³) is the density of the ith snowflake, and N is the total number of snowflakes collected on the plate during the given time interval (Δt_res). Using the average density of the snowflakes during a specific time period and assuming no space between snowflakes or overlapping snowflakes, the new snow accumulation rate $\dot{H}$ (mm h⁻¹) is then

$$\begin{array}{}\text{(18)}& \dot{H}=k{\displaystyle \frac{\mathrm{\Delta }m}{{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}{A}_{\mathrm{hp}}\mathrm{\Delta }{t}_{\mathrm{res}}}}.\end{array}$$

Finally, the total accumulated snow over H (mm) over a given time interval Δt_res (h) is given by $H=\dot{H}\times \mathrm{\Delta }{t}_{\mathrm{res}}$. Note that the bulk density of a fresh snowpack and the height of snowpack can differ from the average density of individual snowflakes (${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}$) and H, respectively, because snowflake settling and compaction on the ground depend on considerations such as their settling characteristics, fall angle, wind speed, the structure of snowflakes, and ambient temperature. We do not account for these processes in the calculation of the volume of freshly fallen snow layers as the impacts are largely unknown. Hence, these variables are proxies for those in the actual snowpack.

3 Experimental methods

Two laboratory experiments and one field experiment were designed to calibrate and validate the MS method for determining snowflake density. The first lab experiment was used to estimate v_melt of ice particles for a given set of environmental conditions and validate the density measurements of ice particles. The second lab experiment investigated the impact of environmental factors on v_melt. A field experiment was conducted at Alta Ski Area's mid-Collins snow-study plot to provide an opportunity to validate the MS method against manual measurements, ultrasonic snow depth sensors, and a weighing gauge using an industry standard method.

3.1 Laboratory experiment method and validation

As illustrated in Fig. 3, laboratory experiments were conducted using a DEID disdrometer, a temperature and relative humidity sensor (Vaisala HMP155, sample rate of 1 Hz, accuracy of ±5 % relative humidity, temperature of ±1 °C), a high-precision gravity scale (Sartorius model ENTRIS64-1S with a readability of 0.1 mg and a standard deviation of 0.1 mg), a micropipette (accuracy of 1 % and maximum accuracy of 1.2 µL at the highest volume), a silicon mold, and a freezer with a minimum temperature of −37 °C. The precise setup of the DEID is modifiable, but for this study, the thermal camera measured surface temperatures of a hotplate with a dimension of ≈ 9 cm × 6 cm at 15 fps with 531 pixels × 362 pixels, which yields a spatial resolution of about 0.2 mm per pixel on the plate. The thermal camera looked vertically downward at the hotplate at an angle. As a result, the maximum error in area measurement is 1.6 %, and the error in the area was corrected using a custom-made function based on the height and angle of the thermal camera, and details are given in Appendix D. Continuous thermal-camera imagery (recorded at 15 Hz) provided all relevant parameters for calculating hydrometeor mass, except for the effective thermal-conduction coefficient between the hotplate and hydrometeor (κ), which was determined through laboratory calibration. We used an Infratec thermal camera that writes out infrared binary (IRB) files that store each pixel's absolute temperature. In post-processing, a grey-scale thermal image ranging in intensity from 0 to 255 is created from the IRB files based on a preset infrared temperature range. To determine Δt_melt and ΔT_melt, the temperature range used was [−40, 0] °C, and for the m, A_c, Δt_evap, and ΔT_evap measurements, the temperature range used was (0, 85] °C, as illustrated in Fig. 3c. In a table-top experiment, the DEID was operated at 85 °C, as determined using the thermal camera. To measure v_melt using h and Δt_melt, the DEID was placed in a 0.25 m per side open-topped cubic enclosure within an environment with a near-zero wind-tunnel wind speed of 0.02 m s⁻¹, a constant ambient temperature of 18 °C, and a constant relative humidity of 38 %. A total of 80 hemispherically shaped ice particles with a range of known masses and volumes were made in a laboratory freezer using a micropipettor by applying a distilled water droplet to a flat silicon mold. A side view of the thermal image of the hemispherical ice particle sample is provided in Fig. 3b. Droplet volumes of 5, 10, 20, 30, 40, 50, 60, and 70 µL, with 10 samples per volume, were used to manufacture the ice particles. The height of the ice particles h was measured in the freezer as illustrated in Fig. 3d. A U-shaped rigid metallic frame was mounted on a translational stage. A millimeter ruler was attached to one arm of the U frame, and a laser pointer was affixed perpendicularly to the second arm facing towards the ruler on the opposite arm in such a manner that when turned on, the pointer's laser beam would strike a spot on the ruler bar. The pointer was horizontally aligned independently using a precision bubble level. The ice particle was thus positioned between the two arms of the U frame. The ruler scale's purpose was to obtain the vertical height of the ice particles. The U frame was vertically displaced when the ice particle occluded the laser beam and did not reach the ruler scale. Upon emerging at the top of the ice particle, when the beam spot hit the ruler scale, the reading was taken as illustrated in Fig. 3d. Furthermore, ice particle mass was also measured with a gravimetric scale prior to its application on the hotplate. The laboratory ice particle sample dimensions are summarized in Table 1.

Figure 3 Details of the DEID MS measurement technique. (a) A hemispherical ice particle applied to the hotplate seen as a bright circular region after melting alongside a rectangular piece of Kapton tape (ϵ ≈ 0.95) used to measure the hotplate surface temperature (T_p). (b) Surface-temperature contour plot of the side view of an ice particle obtained using the thermal camera at t=0 for a temperature range of [−40,0] °C. (c) Schematic of the DEID showing the imaging of melting and evaporating particles, respectively. The black and white contrast of the ice and water particles is optimized by adjusting the camera's temperature range. (d) Schematic of the ice particle height measurement apparatus.

Table 1 The shape of the sample ice particles made in a laboratory: 𝒱_pipette is the volume of water applied by the pipette, A_c is the cross-sectional area of ice particles measured by the thermal camera on the hotplate, and R is the maximum thickness (radius) of ice particles determined from the laser pointer system shown in Fig. 3.

Individual frozen droplets were placed on the hotplate, and the cross-sectional area A_c, m, ΔT_melt, ΔT_evap, Δt_melt, and Δt_evap were measured using the DEID from Eq. (8). v_melt was then calculated using two different formulas $v_{\mathrm{melt}}=h/\mathrm{\Delta }{t}_{\mathrm{melt}}$ and using Eq. (10), where h was determined from the laser pointer system.

3.2 Environmental impacts on v_melt measurement: wind-tunnel experiments

For any given ice particle mass m, independent of the rate of melting or evaporation, the ≈ m(L_f+L_v) constant total quantity of energy is required both to melt and to evaporate a particle from the hotplate. However, the conductive heat rate from the hotplate to the ice particle is a function of environmental conditions through the temperature difference ΔT. To determine how environmental variability affects measurement of the melting velocity v_melt, a portable wind tunnel was set on one side of the DEID's hotplate, allowing different wind velocities to pass over the hotplate and ice particles. A pitot-static probe and an automated weather station measured air speed, ambient temperature, and relative humidity. Ice particles of 60 µL (0.06 g) with radius R=3.25 mm and area $A_{\mathrm{p}}=\mathrm{30.13}\times {\mathrm{10}}^{-\mathrm{6}}$ m² were placed on the DEID's hotplate. Three experiments were performed. First, the hotplate temperature and the ambient relative humidity were maintained constant at 85 °C and 38 %, respectively, while the wind tunnel was adjusted for air speeds of 0.0, 1.5, 3.0, 4.7, 5.9, and 8.3 m s⁻¹. Second, ice particle experiments were performed for surface plate temperatures of 65, 85, and 95 °C for near-zero wind speed and 38 % relative humidity. Finally, the relative humidity was varied to cover 38 %, 68 %, and 91 % with near-zero wind speed and a constant hotplate temperature of 85 °C.

3.3 Density of individual frozen saltwater particles

To test the MS method on a wide range of particle densities, a method was required to produce particles of varying densities. This was done by creating frozen droplets by adding sodium chloride in a distilled water solution (see Table 4 for the percentage of sodium chloride). The densities of these particles were then measured using the DEID and the same methods described above for the pure frozen-water tests.

3.4 Field validation experiments

Data were obtained from field experiments conducted during the winter between October 2020 and April 2021 comprising 17 snowfall events in the upper Little Cottonwood Canyon, Utah, USA, at the Alta Ski Area's mid-Collins snow-study plot (Alcott and Steenburgh, 2010) (40.5763° N, 111.6383° W; 2920 m above sea level). A 10 m crank-up measurement tower at the site included a DEID (sampling rate of 15 Hz) for measurement of microphysical properties of snowflakes and a 3D sonic anemometer (Campbell Scientific, Inc. CSAT3, sampling rate of 20 Hz and accuracy of ±0.05 m s⁻¹). In addition to the DEID, a particle imaging system, consisting of a laser sheet with a sampling volume of 10 cm × 18 cm × 7 cm, was simultaneously deployed and oriented normal to the viewing angle of a Nikon D850 single-lens reflex (SLR) camera as shown in Fig. 4a. The SLR camera recorded 1920 pixel × 1080 pixel images at a spatial resolution of ≈ 160 µm per pixel at 120 fps within a vertical laser sheet created using three 10 W 520 nm diode lasers and a collimator lens. The laser beam spread angle of ≈ 6.8° allowed for a light sheet with near-constant thickness of ≈ 7 cm throughout the region of interest. A single-focal-length Nikon AF-S VR Micro-Nikkor 105 mm $f/\mathrm{2.8}$ G IF-ED lens permitted a depth of field greater than the thickness of the laser light sheet. The DEID was deployed 2 cm below the lower end of the laser sheet, which permits measurement of the microphysical properties of snowflakes that pass the laser sheet and fall on the hotplate. A Vaisala HMP155 temperature and relative humidity sensor (1 Hz sampling rate) was also located on the tower and maintained at a height of approximately 1.5 m above the new snow level. At approximately the same height, a Campbell Scientific, Inc. CSAT3 3D sonic anemometer was deployed (sampling rate of 20 Hz).

Figure 4 (a) Field experiment setup for measurements of microphysical properties of snowflakes and snowflake visualization. The experimental setup consists of (1) a 20 m tower, (2) a thermal camera, (3) a hotplate, (4) three 10 W lasers and an optical lens, (5) a D850 Nikon SLR camera, (6) a 3D sonic anemometer, (7) a data logger and computer, and (8) a relative humidity and temperature sensor. (b) Top-view schematic showing the co-location of the laser system and DEID hotplate.

Images from the SLR camera were combined with mass measurements from the DEID to compute snowflake density in the field and validate the MS method. With this SLR–DEID method, the geometrical volume of each free-falling snowflake was estimated using images from the SLR camera. The mass of each hydrometeor was determined by following individual snowflakes through the laser sheet until they hit the DEID hotplate. This method was applied to approximately 1000 snowflakes. Selected thermal images of aggregate snowflakes and graupel on the hotplate are illustrated in Appendix B.

3.5 Hydrometeor density calculation exploiting concurrent imagery during their fall

The geometrical volume 𝒱_SLR of a snowflake can be estimated independently of the DEID by imaging falling snowflakes using the particle tracking system discussed above. In this case, the density is determined from

$$\begin{array}{}\text{(19)}& {\mathit{\rho }}_{\text{SLR-DEID}}={\displaystyle \frac{m}{{\mathcal{V}}_{\mathrm{SLR}}}},\end{array}$$

where snowflake mass m is determined with the DEID from Eq. (8).

We categorized five snowflake habits based on the international classification for seasonal snow on the ground (Fierz et al., 2009; Praz et al., 2017): planar crystal (combining stellar and plates), graupel (combining hail and graupel), columnar crystal, aggregate (combining irregular crystal), and small particles. Graupel and small-particle crystals were classified based on size. Each snowflake category contained approximately 200 samples. Taking advantage of how snowflakes rotate while falling, a single camera with multiple images was found to represent a 3D picture of a snowflake and provide geometrical volume more accurately than using a single image (Li et al., 2022). Five sequentially selected images of each crystal type are illustrated in Fig. 5a. A schematic showing how snowflake volume is computed is illustrated in Fig. 5b and formulated in Table 2. The volume of planar crystals was approximated as a disk; hence, ${\mathcal{V}}_{\mathrm{SLR}}=A_{max}{D}_{min}$. Here, A_max is the maximum area of a planar crystal in a 2D plane among all the images that are visible to the camera, and D_min is the minimum dimension (representing the thickness of the disk) in the 2D plane as indicated in Fig. 5a, column I, row V. The volume of graupel was estimated as a sphere, ${\mathcal{V}}_{\mathrm{SLR}}=\mathit{\pi }/\mathrm{6}{D}_{\mathrm{Reff}}^{\mathrm{3}}$, where D_Reff is the effective circular diameter estimated using the SLR camera imaging the snowflake falling in the air and estimated as $D_{\mathrm{Reff}}=\sqrt{\frac{\mathrm{4}}{\mathit{\pi }}\stackrel{\mathrm{‾}}{A}}$. $\stackrel{\mathrm{‾}}{A}$ is the 2D mean area of all images as illustrated in Fig. 5a, column II. The volume of aggregates was estimated by fitting an ellipsoid such that ${\mathcal{V}}_{\mathrm{SLR}}=\frac{\mathit{\pi }}{\mathrm{6}}{D}_{max}{D}_{min}{D}_{\mathrm{v}}$, where D_max, D_min, and D_v are the lengths of three mutually perpendicular ellipsoid axes as illustrated in Fig. 5b. The volume of columnar crystals (solid cylinder) was estimated from ${\mathcal{V}}_{\mathrm{SLR}}=\frac{\mathit{\pi }}{\mathrm{4}}{\stackrel{\mathrm{‾}}{D}}_{min}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{D}}_{max}$, where ${\stackrel{\mathrm{‾}}{D}}_{max}$ and ${\stackrel{\mathrm{‾}}{D}}_{min}$ are the average of the maximum and minimum dimensions of five images as illustrated in Fig. 5a, column III. The volumes of small particles were estimated using a spherical volume ${\mathcal{V}}_{\mathrm{SLR}}=\frac{\mathit{\pi }}{\mathrm{6}}{\stackrel{\mathrm{‾}}{D}}_{max}^{\mathrm{3}}$, where ${\stackrel{\mathrm{‾}}{D}}_{max}$ is the average of the maximum dimension of five small-particle images (Fig. 5a, column V).

Figure 5 (a) For five snowflake types, five images of each snowflake separated due to rotation just prior to settling on the DEID's hotplate. (b) Volume measurement method.

Table 2 Classification of snowflakes and corresponding volume estimation based on crystal geometry.

3.6 Manual measurements of the snowpack: bulk density snowpack calculations

The mean bulk density of a fresh snowpack (${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{m}}$) can also be determined using manual measurements of the ratio of the snow water equivalent depth (SWE) to the new snow depth (H). This can be done with the DEID by recalling that $\mathrm{SWE}=k\mathrm{\Delta }m/{\mathit{\rho }}_{\mathrm{w}}{A}_{\mathrm{hp}}$ and $H=k\mathrm{\Delta }m/{\mathit{\rho }}_{\mathrm{s}}^{\mathrm{m}}{A}_{\mathrm{hp}}$, where

$$\begin{array}{}\text{(20)}& {\mathit{\rho }}_{\mathrm{s}}^{\mathrm{m}}={\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\mathrm{SWE}}{H}}.\end{array}$$

Here, ρ_w is the density of water of 1000 kg m⁻³. At the Alta-Collins snow-study plot these measurements were obtained every 12 h. Note that changes within the snowpack due to processes such as densification, heat transfer, wind shear, etc. are not considered here as our analyses are limited to the consideration of freshly fallen snow.

For further comparison, the average snowpack density can also be estimated at hourly intervals based on measurements obtained from the ETI Instrument Systems Noah II precipitation weighing gauge sensor (SWE_ETI) and the snow depth from the Campbell Scientific, Inc. SR50 ultrasonic snow depth sensor. The ETI and SR50 sensors were deployed 4 m from the DEID at the Alta-Collins site. A windshield was deployed around the ETI bucket to increase catchment efficiency. The ETI reported SWE measurements once every hour with a resolution, threshold, and accuracy of 0.25, 0.25, and ±0.25 mm, respectively. The SR50 sensor recorded snow depth every hour to provide running totals of snow depth. The measurement range of the ultrasonic sensor was 0.5 to 10 m with an accuracy of 0.4 % and a resolution of 0.1 mm. Raw DEID data sampled at a rate of 15 Hz were integrated to produce hourly measurements for comparison with ETI data. ETI and SR50 data were collected throughout the winter of 2020, but data from 07:00 UTC on 12 December to 19:00 UTC on 12 December 2020 were used to compare SWE and snow accumulation with the DEID obtained using the MS density measurement method.

4 Results

4.1 MS density method laboratory calibration

In order to use the MS method for determining density, the melting calibration constant, c_melt, and the calibration constant, c, in Eqs. (10) and (12), respectively, must first be determined empirically. To do this, 80 ice particles with different sizes and masses were applied to the DEID hotplate. Furthermore, the experiments were conducted in an environmentally controlled chamber with the air temperature fixed at 18 °C, near-zero wind velocity (0.05 m s⁻¹), a hotplate temperature of 85 °C, and 38 % relative humidity. Results from experiments conducted at these standard conditions are identified with subscript 0. The variables Δt_melt, Δt_evap, ΔT_melt, and ΔT_evap for each particle were determined using the thermal camera, while h was measured directly using the laser pointer system. The measured values of h and Δt_melt are plotted in Fig. 6a. The slope of the h–Δt_melt curve is v_melt, which is approximately constant ($v_{\mathrm{0},\mathrm{melt}}=\mathrm{2.11}\pm \mathrm{0.10}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹) because the experiments were performed with ice particles in an environmentally controlled chamber where the average measured value of ΔT_0,melt was also found to be approximately constant (3.15 ± 0.02 K)​​​​​​​​​​​​​​.

Figure 6 (a) Ice particle height as a function of Δt_melt. The slope of this line determines v_0,melt, the melting speed under fixed standard conditions. That is, ice particles of different maximum thicknesses [0.22, 5.6] mm and their melting time at a plate temperature of 85 °C, near-zero wind velocity (0.05 m s⁻¹), and 38 % relative humidity. (b) Plot of v_melt versus c_meltΔT_melt illustrating the validity of Eq. (10). v_melt is determined directly from measurements of particle maximum thickness and melting time, Δt_melt, and then compared to c_meltΔT_melt for a range of environmental conditions.

The measured v_melt and ΔT_melt for each particle can be substituted into Eq. (10) to solve for c_melt, and then the constant c can easily be determined from Eq. (13). This was done, and the results were averaged over all 80 samples, yielding $c_{\mathrm{melt}}=\mathrm{6.69}\pm \mathrm{0.11}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ K⁻¹​​​​​​​ and $c=\mathrm{1.01}\pm \mathrm{0.33}\times {\mathrm{10}}^{\mathrm{5}}$ K s m⁻¹​​​​​​​​​​​​​​. With a derived value of c and particle mass measured with the DEID, the particle density can now be inferred from Eq. (12). This yields an average ice particle density of ${\mathit{\rho }}_{\mathrm{MS}}^{\mathrm{ice}}=\mathrm{919}\pm \mathrm{65}$ kg m⁻³. This is very close to the expected value of ice density at temperatures near 0 °C (i.e., 917 kg m⁻³).

A summary of the following measured parameters for laboratory-created ice particles is presented in Table 3: v_melt, ${\mathit{\rho }}_{\mathrm{MS}}^{\mathrm{ice}}$, m_DEID, m_gravity, and 𝒱_pipette. Here, 𝒱_pipette is the volume of a water droplet, not the ice particle volume. Since these experiments were conducted in an enclosure where environmental variability was negligible, the effect of convective cooling on the measurement of m, A_c, Δt_melt Δt_evap, and ΔT_evap did not play a role. However, in nature, winds can affect Δt_melt Δt_evap, ΔT_evap, and v_melt.

Table 3 Mean and standard deviation of the mass of 80 experimentally manufactured ice particles, with the water droplet volume applied by pipette (𝒱_pipette) prior to freezing, determined using a gravimetric scale (m_gravity) and the DEID (m_DEID). The density of ice particles estimated using the DEID MS method – ${\mathit{\rho }}_{\mathrm{MS}}^{\mathrm{ice}}$. The melting velocity under the standard conditions described in the text – v_0,melt (m s⁻¹).

To determine how environmental variability affects v_melt and to generalize the validation of Eq. (10), 60 µL (0.06 g) ice particles with thickness R=3.25 mm and area $A_{\mathrm{c}}=\mathrm{30.13}\times {\mathrm{10}}^{-\mathrm{6}}$ m² were placed on the DEID's hotplate, and the wind speed was varied from 0.0 to 8.3 m s⁻¹, relative humidity varied from 38 % to 91 %, and plate temperature varied from 65 to 95 °C. v_melt was computed using direct measurements as $v_{\mathrm{melt}}=R/\mathrm{\Delta }{t}_{\mathrm{melt}}$. It was also computed using Eq. (10; i.e., v_melt=c_meltΔT_melt. Results are shown in Fig. 6b. The coefficient of determination R² between v_melt computed using the two methods is 0.99, and the RMS error is 3.46 × 10⁻⁶ m s⁻¹. Using Eq. (12), the average measured density of 80 ice particles with different shapes and sizes, formed by pipetting water droplets onto arbitrary surfaces with different contact angles under the above-listed environmental conditions and plate temperatures, was 928 ± 56 kg m⁻³, which shows the MS method works for a wide range of environmental conditions within experimental uncertainty.

4.2 Density of individual frozen saltwater particles

To test the MS method on a wider range of particle densities, frozen saltwater particles with different salt (NaCl) concentrations were applied to the DEID hotplate. The estimated density of the frozen saltwater particles calculated using the MS method agrees with the density expected from the percentage of salt in the solution (ρ_theoretical) to within 3 %. The results are summarized in Table 4.

Table 4 The density of frozen saltwater particles measured using the MS method and determined theoretically based on a saltwater ratio by weight (percentage of salt and water).

4.3 Field evaluation of the MS method for snow particles of different types

We collected data at the Alta Ski Area's mid-Collins snow-study plot from 7 November 2020 to 27 April 2021. During this time, a snow accumulation of 12.35 m and an SWE accumulation of 1.38 m were observed with the DEID. The ambient air temperature varied from −21 to 2 °C, relative humidity varied from 64 % to 97 %, and wind speed varied from 0.2 to 8 m s⁻¹. Generally, the observed densities of freshly fallen individual snowflakes varied from 9 to 495 kg m⁻³. The average densities of each storm varied from 35 to 115 kg m⁻³. Figure 7 shows estimated snowflake densities using both the SLR–DEID method and the MS method for five types of snowflakes. A comparison between the MS and SLR–DEID density methods for five crystal types is summarized in Table 5. The coefficient of determination between the two methods is the highest for small particles and graupel and the lowest for aggregates. The measured size (D_eff) of each type of crystal is summarized in Table 5. The uncertainty that arises between the two methods for aggregate snowflakes may be due to the SLR–DEID method used to estimate the geometrical volume because aggregates have a more irregular shape than the other snow crystal types. For all snowflakes, the mean estimated density is 131 ± 83 kg m⁻³ using the SLR–DEID method and 142 ± 87 kg m⁻³ using the MS method, yielding an uncertainty of 3.9 % and an R² of 0.95.

Figure 7 Measured density of a range of snowflake types using the particle imaging system compared with values obtained using the DEID MS method with associated coefficients of determination, slopes, and RMSE errors.

Table 5 Comparison between two density methods, MS and SLR–DEID, of five types of snow crystals. Range of D_eff for each type of crystal.

4.4 Frequency distribution of individual hydrometeor densities

Figure 8a shows a probability distribution function (PDF) for the densities of individual snowflakes using the MS method applied to data acquired at the Alta-Collins snow-study plot for seven snow storms selected to encompass a broad range of environmental conditions: mean ambient temperatures [−13.45, −4.82] °C, mean wind speeds [0.30, 0.89] m s⁻¹, and mean relative humidities [72 %, 91 %] as listed in Table 6. The kurtosis (Kr) and skewness (Sk) of the normalized density distribution functions vary from 2.02 to 2.42 and from 0.01 to 0.10, respectively. We found that the PDFs of snow density are symmetric with respect to the mean (Sk=0.01) and near Gaussian (Kr=2.41) when the ambient temperature is the lowest, while Sk=0.10 and Kr=2.02 when the ambient temperature is the highest, which is shown in Table 6. Figure 8b includes results from snowflake densities computed assuming a spherical-particle volume but also using the DEID, as done in Rees et al. (2021). The spherical assumption underestimates snowflake density by a factor of ≈ 1.5 compared to the MS method.

Figure 8 (a) Probability distribution functions of the density of individual snowflakes measured using the MS method for seven storms. (b) Comparison between the densities determined with the two methods employing DEID measurements alone, the MS method, and the spherical-particle method. The solid line, dashed line, and dotted line are the mean densities of a storm using a spherical method, an MS method and manual measurements, respectively. Note that the MS method and manual measurement are compared without knowing the information about snowpack, as mentioned earlier.

Table 6 Summary of DEID-derived snowflake parameters for a series of storms in Alta, Utah. Mean diameter ${\stackrel{\mathrm{‾}}{D}}_{\mathrm{eff}}$, mean density ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}$, and total number of snowflakes N captured during seven storms with a mean ambient temperature ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{amb}}$, mean wind speed $\stackrel{\mathrm{‾}}{U}$, mean relative humidity $\stackrel{\mathrm{‾}}{RH}$, total snow depth, and total accumulated snow water equivalent for seven storms using the DEID. Kr is the kurtosis and Sk the skewness of the density distribution.

4.5 Validation of SWE measurements

SWE determined with the DEID can be compared to manual measurements collected at the Alta-Collins study plot. Since manual measurements are made infrequently at intervals of 12 h, the comparisons are made on a storm-by-storm basis as shown in Fig. 9. The relationship between DEID observations and the bulk standard manual measurement techniques is shown in Fig. 9a. A best-fit relationship between the two methods yields an R² of 0.994 with a slope of 0.94 ± 0.04.

Figure 9 (a) SWE accumulation from DEID and manual measurements for 17 storms. Each data point represents an individual storm. Manual measurements were made every 12 h (11:00 and 23:00 UTC), and the DEID was sampled at 15 Hz. (b) Time series of SWE accumulation and SWE rate measured using the DEID and ETI gauge. Each DEID data point represents a 1 min average and each ETI gauge data point a 1 h average.

The accumulated SWE integrated over 1 min intervals is compared to ETI data in Fig. 9b. DEID SWE accumulation observations match those from the ETI gauge to within ±6 % over the 12 h measurement period (storm duration). SWE accumulation measured by the DEID is slightly higher than that obtained by the ETI because the minimum resolution of the ETI is 0.254 mm, whereas the minimum DEID resolution is 0.001 mm (Singh et al., 2021). Furthermore, the ETI gauge has been shown to undercatch snowflakes compared to the DEID under high-wind conditions (Singh et al., 2021).

4.6 Validation of snow depth measurements

Snow accumulation (H) was computed with the DEID using Eq. (18) and $H=\dot{H}\times \mathrm{\Delta }{t}_{\mathrm{res}}$, where 1 min average MS density was used. The total snow accumulated in each storm measured using the DEID MS density method compares well with manual measurements obtained every 12 h at the Alta-Collins snow-study plot with an R² value of 0.983 and a slope of 1.12 ± 0.07 as shown in Fig. 10a. A 12 h period with snowfall between 07:00 and 19:00 UTC on 12 December 2020 is shown in Fig. 10b. There is also good agreement to within ±5 % between snow accumulation measurements obtained from the DEID using the MS density method and those obtained using an ultrasonic snow depth sensor obtained once per hour. The bulk density of a fresh snowpack can differ from the average density of individual snowflakes and the 1 min average because snowflake settling and compaction on the ground depend on considerations such as their settling characteristics, fall angle, wind speed, the structure of snowflakes, and the ambient temperature. We do not account for these processes in the calculation of the volume of freshly fallen snow layers as the impacts are largely unknown. Nonetheless, the bulk density of the snowpack measured during 17 storms using the DEID MS method can also be compared with manual gravimetric snow density measurements of SWE depth (SWE_m) and snow depth (H) (Eq. 20). The R² value relating the DEID and manual bulk density measurements is 0.88 with a slope of 0.90 ± 0.15, as shown in Fig. 11. The implication of these two comparisons, somewhat surprisingly, is that the DEID reproduces measurements of freshly fallen snowpack density and accumulation, made with more traditional techniques, without considering the quite complex physics of how individual snowflakes pack and stack.

Figure 10 (a) For 17 storms, comparison between snow accumulation from the DEID obtained using the MS density method and from manual measurements made every 12 h (11:00 and 23:00 UTC). Each data point represents a storm. (b) Time series of snow accumulation and snowfall rate measured using 1 min averaged DEID results and 1 h averaged ultrasonic snow depth sensor measurements during a storm on 12 December 2020.

Figure 11 For 17 storms, comparison between bulk density of snowpack from the DEID obtained using the MS density method and from manual gravimetric snow density measurements made every 12 h (11:00 and 23:00 UTC). Each data point represents a storm.

5 Conclusions

Automated determination of the density of individual snowflakes has been a long-standing challenge. In this study, we present a novel method for accomplishing this goal that exploits a new hotplate instrument, the Differential Emissivity Imaging Disdrometer or DEID, which has previously been shown to be capable of obtaining highly accurate direct measurements of particle mass. Particle-by-particle density estimates are obtained from measurements of particle mass, particle contact area onto the hotplate, and an estimate of the particle's effective thickness using melting speed and melting time. A particle's effective thickness normal to the hotplate is a product of the concept of a melting speed and melting time from which, in combination with a particle's contact area, an estimate of particle volume and hence its density can be obtained. For individual hydrometeors, this melting speed method was validated using laboratory ice particles of known density showing a maximum uncertainty of 6.3 %, as well as with videos from field measurements of a range of naturally falling snowflakes with an uncertainty of 3.7 %. DEID observations at a high-elevation mountain site of snow water equivalent (SWE) accumulation, snow depth accumulation (H), and bulk snow density (${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{MS}}$) from 17 storms taken at the Alta-Collins snow-study plot at the Alta Ski Area during the winter of 2020–2021 agreed well with traditional manual technique measurements with R² values of 0.994, 0.983, and 0.88 and slopes of 0.94 ± 0.04, 1.12 ± 0.07, and 0.90 ± 0.15, respectively, independent of environmental conditions, including wind speed, ambient temperature, relative humidity, and hotplate temperature. We acknowledge that these bulk snow results are surprisingly good given the complexity of the processes associated with snow once it hits the ground. We speculate that the snow depth accumulation and bulk snow density results imply minimal snowflake overlap in the snowpack. That is, the sum of the heights of the snowflakes is equal to the total depth. It also implies that the fresh snowpack minimally settles down during the short duration (12 h) associated with the manual measurements. In reality, there is overlap, but on average, the snowflakes are approximately non-overlapping. Our results imply that the density of the snowpack is equal to the average density of individual snowflakes, which means that the snowflakes “pack” with a density similar to the density of the snowflake. We note that sometimes this is not true, for example, when we have an error of 15 %.

Using the melting speed method, the utility of the DEID design can be extended from hydrometeor mass measurement to measurement of the density of irregularly shaped hydrometeors, in real time with high accuracy. Such information, whether taken on a particle-by-particle basis or assessed as a cumulative bulk quantity, can be applied to high-resolution measurement of vertical density variability in the snowpack (e.g., Morrison et al., 2023) – critical for the assessment of its stability – and to studies of snow metamorphism, the assessment of transitions between rain and snow, determinations of rates of hydrometeor settling and its response to turbulence, and the scattering of visible light and radar pulses by snow particles.

Appendix A: Latent heat of vaporization calculation

The latent heat of vaporization of water depends on temperature and may be written as $L_{\mathrm{v}}\left(T\right)=(\mathrm{2.501}-\mathrm{0.00237}\times T)\times$10⁶ J kg⁻¹, where T is in degrees Celsius (Stull, 2012). In our case, water droplets evaporate at temperatures from 0 to 85 °C. A sample time series of the temperature of a water droplet during evaporation is plotted in Fig. A1. The estimated L_v based on the time series of temperature of water droplets is $\mathrm{2.32}\pm \mathrm{0.02}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹.

Figure A1 Sample time series of the temperature of the melted portion of an ice particle after being placed on the hotplate. A 60 µL water droplet volume was used to make the ice particle shown here.

Appendix B: Relation between the melting and evaporation statistics

Using Eq. (7), the mass associated with the ijth pixel can be estimated separately during melting and evaporation as shown in Eqs. (B1) and (B2) below, respectively.

$$\begin{array}{}\text{(B1)}& {\displaystyle }{m}_{ij}={\displaystyle \frac{\mathit{\kappa }}{C_{\mathrm{ice}}{T}_{\mathrm{0}}+L_{\mathrm{f}}}}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{ij,\mathrm{melt}}}{\int }}\mathrm{\Delta }{T}_{ij,\mathrm{melt}}\left(t\right)A_{ij}\left(t\right)\mathrm{d}t\text{(B2)}& {\displaystyle }{m}_{ij}={\displaystyle \frac{\mathit{\kappa }}{C_{ij,\mathrm{liq}}{T}_{\mathrm{p}}+L_{\mathrm{v}}}}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }{t}_{ij,\mathrm{evap}}}{\int }}\mathrm{\Delta }{T}_{ij,\mathrm{evap}}\left(t\right)A_{ij}\left(t\right)\mathrm{d}t\end{array}$$

Here, A_ij is the contact area of the single pixel, which is constant. By using Eqs. (B1) and (B2), we may write

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{\Delta }{T}_{ij,\mathrm{melt}}\mathrm{\Delta }{t}_{ij,\mathrm{melt}}}{L_{\mathrm{ff}}}}={\displaystyle \frac{\mathrm{\Delta }{T}_{ij,\mathrm{evap}}\mathrm{\Delta }{t}_{ij,\mathrm{evap}}}{L_{\mathrm{vv}}}}.\end{array}$$

In Eq. (B3), $L_{\mathrm{ff}}=C_{\mathrm{ice}}{T}_{\mathrm{0}}+L_{\mathrm{f}}$ and $L_{\mathrm{vv}}=C_{\mathrm{liq}}{T}_{\mathrm{p}}+L_{\mathrm{v}}$.

Equation (B3) can be re-written for all pixels associated with hydrometers using double summation.

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}& \sum _i\sum _j{\displaystyle \frac{\mathrm{\Delta }{T}_{ij,\mathrm{melt}}\mathrm{\Delta }{t}_{ij,\mathrm{melt}}}{L_{\mathrm{ff}}}}\\ & =\sum _i\sum _j{\displaystyle \frac{\mathrm{\Delta }{T}_{ij,\mathrm{evap}}\mathrm{\Delta }{t}_{ij,\mathrm{evap}}}{L_{\mathrm{vv}}}}\end{array}\end{array}$$

Based on experimental tests, we assume $\mathrm{\Delta }{t}_{ij,\mathrm{melt}}\approx \mathrm{\Delta }{t}_{\mathrm{melt}}$ and $\mathrm{\Delta }{t}_{ij,\mathrm{evap}}\approx \mathrm{\Delta }{t}_{\mathrm{evap}}$ (Eq. B4), which yields

$$\begin{array}{}\text{(B5)}& {\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{melt}}\mathrm{\Delta }{t}_{\mathrm{melt}}}{L_{\mathrm{ff}}}}\approx {\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{evap}}\mathrm{\Delta }{t}_{\mathrm{evap}}}{L_{\mathrm{vv}}}}.\end{array}$$

Experimentally determined values for ΔT_melt, Δt_melt, ΔT_evap, and Δt_evap for different ice particle masses and sizes for a large range of environmental conditions are shown in Fig. B1. As expected, the slope between the two terms agrees to within 6.5 % of the ratio of L_vv to L_ff. With the DEID, Δt_melt and Δt_evap are directly measured using a thermal camera and counting the number of frames between the first frame when an ice particle hits the hotplate and last frame when the particle has completely melted or evaporated.

Figure B1 Experimentally measured Δt_melt, ΔT_melt, Δt_evap, and ΔT_evap for laboratory ice particles. The slope of the line is quite close to $L_{\mathrm{vv}}/L_{\mathrm{ff}}$. An error bar is associated with 10 samples in each ice particle size.

Appendix C: DEID thermal imagery

In this section, we present example thermal images of different types of snowflakes after melting and during evaporation on the DEID hotplate. Figures C1 and C2 show aggregate and graupel snow particles, respectively. These data were taken at the Alta-Collins snow-study plot on 22 December 2020.

Figure C1 (a) Black and white binary thermal images of an aggregate type of snowflake at different stages of melting and evaporation on the DEID hotplate observed in Alta, UT. (b) Cropped aggregate snowflake images on the hotplate just after melting.

Figure C2 Black and white binary thermal images of a graupel type of snowflake during different stages of melting and evaporation on the DEID heated plate observed at Alta.

Appendix D: Systematic and random error analysis

Direct measurements made by the DEID include snowflake area, temperature, and evaporation time, for which the respective uncertainties are 1.4 %, 0.3 %, and 1.6 %. Due to the inclination of the thermal camera, the maximum error in the measurement of area is 1.6 %, and it is corrected using a custom-made MATLAB function based on the height and angle of the thermal camera. Both x and y direction pixels were corrected using algorithm h_camtan (θ_y+jdθ). Here, h_cam is the height of the camera and ${\mathit{\theta }}_y=\mathrm{90}-{\mathit{\theta }}_H-{\mathit{\theta }}_V/\mathrm{2}$. The angles θ_H and θ_V are the thermal camera's horizontal and vertical angles, $d\mathit{\theta }={\mathit{\theta }}_V/j$, where j is the number of pixels in the y direction. A similar method was used for the x direction. The uncertainties in derived quantities (using a standard propagation of uncertainty analysis), such as mass (m), height (h), and density (ρ_s), are 4.3 %, 2.9 %, and 8.6 %, respectively.

Appendix E: Separate quantification of area and temperature of ice and liquid

During the post-processing of the data, the temperature ranges of the thermal images were set such that only one phase, either ice or liquid, could be seen on the hotplate. Figure E1a shows a grey thermal image where the temperature of some portions is less than or equal to 0, and some are greater than 0. When using the temperature range [−40, 0] °C, only the ice portion is visible, as seen in Fig. E1b. The sum of all visible areas is A_ice(t), and the spatial mean temperature over all those pixels is called T_ice(t). Similarly, when the temperature range (0, 85] °C is used, only the liquid portion is visible, as seen in Fig. E1c. The sum of all visible areas is A_liq(t), and the spatial mean temperature over all those pixels is called T_liq(t).

Figure E1 (a) Schematic illustrating the measurement area (contact area) and temperature of ice and liquid separately. (a) Thermal image of the frozen hydrometeor on the hotplate. The temperatures less than or greater than 0 show the frozen hydrometeor's unmelted and melted portions, respectively. The temperature recorded by the thermal camera for the hotplate is around 0 due to the low emissivity of the plate, where the actual temperature is 85 °C. (b) During the post-processing of the data, the temperature range was set to [−40,0] °C, which allows one to see only the area of the unmelted portion (temperatures less than or equal to 0) of the frozen hydrometeor, and the area of the melted portion (temperatures greater than 0) is 0. T_ice(t) is the mean temperature of all pixels with a temperature less than or equal to 0, and those pixels contrast well with the background, which allows them to be counted easily. (c) During the post-processing of the data, the temperature range was set to (0, 85] °C, allowing us to see only the area of the melted portion (temperature greater than 0) of the hydrometeor. T_liq(t) is the mean temperature of all pixels with a temperature greater than 0, and those pixels contrast well with the background, which allows them to be counted easily.

Appendix F: Image of water droplets on molds

This section presents example images of water droplets on different molds while making laboratory ice particles. Figure F1a and b show water droplets on a flat silicon mold and random surface, respectively. The contact angle between water and the silicon mold is about 90°, and the water droplet looks near hemispherical.

Figure F1 (a) A photo of a water droplet on the silicon mold, taken from the side with the SLR camera. (b) An image of a water droplet on the random surface. Note the difference in contact angle.