The Precipitation Imaging Package: Assessment of Microphysical and Bulk Characteristics of Snow

Abstract

Remote-sensing observations are needed to estimate the regional and global impacts of snow. However, to retrieve accurate estimates of snow mass and rate, these observations require augmentation through additional information and assumptions about hydrometeor properties. The Precipitation Imaging Package (PIP) provides information about precipitation characteristics and can be utilized to improve estimates of snowfall rate and accumulation. Here, the goal is to demonstrate the quality and utility of two higher-order PIP-derived products: liquid water equivalent snow rate and an approximation of volume-weighted density called equivalent density. Accuracy of the PIP snow rate and equivalent density is obtained through intercomparison with established retrieval methods and through evaluation with colocated ground-based observations. The results confirm the ability of the PIP-derived products to quantify properties of snow rate and equivalent density, and demonstrate that the PIP produces physically realistic snow characteristics. When compared to the National Weather Service (NWS) snow field measurements of six-hourly accumulation, the PIP-derived accumulations were biased only +2.48% higher. Additionally, this work illustrates fundamentally different microphysical and bulk features of low and high snow-to-liquid ratio events, through assessment of observed particle size distributions, retrieved mass coefficients, and bulk properties. Importantly, this research establishes the role that PIP observations and higher-order products can serve for constraining microphysical assumptions in ground-based and spaceborne remotely sensed snowfall retrievals.

1. Introduction

Remote-sensing observations, principally radar, and passive microwave measurements, provide essential information for estimating and forecasting the impacts of precipitation, and relevant observing platforms provide coverage both nationally and globally [1,2,3]. However, these observations contain incomplete information about the precipitation, therefore assumptions about hydrometeor properties are required in order to retrieve accurate estimates of mass and rate. This is especially true for snow, as more complicated physical properties drive estimation uncertainties [4,5,6]. These uncertainties primarily result from variability in microphysical characteristics of the populations of particles, i.e., particle size distributions (PSDs) [7,8]. Further uncertainties arise due to variations in individual particle properties and how they are fundamentally related to the remote-sensing detector response [9,10,11,12]. A video disdrometer instrument, such as the Precipitation Imaging Package (PIP), provides information about hydrometer characteristics (e.g., PSDs, particle fall speed, and higher-order-derived products, such as density and rate), and can be utilized to reduce uncertainty and greatly improve estimates of snowfall rate and accumulation. For example, these in situ observations can be applied to constrain microphysical snow particle modeling and through forward-modeling address the challenges inherent to remote-sensing retrievals [13,14].

Investigations of PIP-derived products, coupled with colocated remote-sensing instrumentation and surface meteorology measurements, are active areas of research. Specifically, Tiira et al. [15] employed the PIP to examine connections between snow density and PSDs and fall velocities to better constrain relationships of mass and particle dimensions. Additionally, studies have used observations from the PIP (and its predecessor, the Snowflake Video Imager [16]) to produce mass-dimensional relationships of snow particles to assess the quality of radar reflectivity-derived snowfall rate relationships [13,17,18]. PIP-derived products of snow properties in conjunction with meteorological measurements and reanalysis products have been used to illuminate microphysical characteristics associated with specific types of snowfall events [19,20], and to generate parameterizations to improve space-based retrievals of precipitating snow properties [8,21].

The ability of spaceborne remote-sensing instruments to detect surface snowfall events and provide accurate snow rate retrievals is complicated by various factors, including inherent instrument limitations, microphysical variability, environmental conditions, and surface type [22,23,24]. Seasonal ground-based and spaceborne radar snow estimate comparisons are reasonable when averaged over relatively large spatial domains [4]. However, other studies have indicated 30–40% spaceborne radar overestimates compared to ground measurements, with results strongly geographically and latitudinally dependent [9]. Additionally, the presence of orography adversely impacts both ground and space-based radars, resulting in underestimation biases >70% [25,26]. In situ ground-based observations of snowfall accumulation also suffer from relatively large uncertainties, as intercomparisons have shown that gauges are biased low, with as much as 50% of the accumulation missed when unshielded and 12–27% with a single perimeter fence [27,28,29].

In this work, we demonstrate the quality and utility of two higher-order PIP-derived products: liquid water equivalent snow rate and an approximate measure of volume-weighted density called equivalent density. These derived products demonstrate the ability of the PIP to obtain accurate estimates of snow accumulation when compared to proven ground-based observation methods and established retrievals of mass. The PIP snow rate and equivalent density retrievals are assessed through intercomparison with established retrieval methods by von Lerber et al. [17] and Wood et al. [30] and through evaluation with colocated ground-based observations. We achieve quality intercomparisons of snow events by leveraging the unique and comprehensive dataset collected at the National Weather Service Marquette, Michigan office [19]. The presented work separates snow events into categories of low or high snow-to-liquid ratio (SLR) and will examine the differences using observations from the PIP and bulk snow characteristics from the PIP higher-order-derived products and the von Lerber et al. [17] and Wood et al. [30] retrievals. We will first describe the data used in this work, as well as the methods employed to obtain estimates of mass and snow rates (Section 2). PIP-derived data products will also be described, with a focus on the empirically derived equivalent density parameterization, as it has not before been published (Section 2.2.1). Section 3 will outline and discuss intercomparison results for the low and high SLR snow events, starting with microphysical and mass characteristics, then focusing on the snow rate and density results.

2. Data and Methods

2.1. Location and Instrumentation

This study utilizes an enhanced precipitation instrument suite located at the Marquette, Michigan (MQT) National Weather Service (NWS) office [19]. This long-duration dataset consists of colocated remotely sensed and in situ observations of precipitation from January 2014 through the present. This region receives approximately 250 to 500 cm of measured snow accumulation per year, half of which is produced by Lake-effect Snow (LeS) processes [19]. Because this area receives both synoptically and surface-forced (e.g., LeS and orographic) snowfall, there is a wide range of observed snow event environmental conditions and snow-to-liquid ratios (SLRs).

2.1.1. Precipitation Imaging Package

The Precipitation Imaging Package (PIP) is a NASA-developed video disdrometer, which evolved from the Snowflake Video Imager (SVI) [16]. The PIP is inexpensive and easily deployable and plays a major role in the NASA Global Precipitation Measurement (GPM) mission ground validation (GV) campaigns [3]. The instrument consists of a high-speed video camera (380 frames per second) with a charge-coupled device image sensor. This camera is aimed at a bright halogen lamp at 2 m distance, with a focal plane at 1.3 m distance. The PIP has a field of view of 64 by 48 mm, an image resolution of 0.1 by 0.1 mm, and a depth of field dependent on the particle diameter. The PIP setup is unique in that precipitation particles are relatively unimpeded by the instrument itself, as the design provides an open volume between the camera and the lamp [16]. This is particularly important for accurately observing and quantifying snowfall, as measurements from other widely used precipitation gauges, i.e., Pluvio, Geonor, etc., will be adversely impacted by even moderate winds [27,28,29,30,31].

The PIP camera records videos of the hydrometeor shadows as they fall through the observation volume in front of the light source. Image processing software then identifies each particle individually and is able to discern associated characteristics from these observations (see Newman et al. [16] for details). For example, the PIP records shape characteristics, such as the particle equivalent area diameter, which is the diameter of a circle for which the area equals that of the particle projected upon the plane of view. Additionally, for particles that are observed in more than one frame, both the horizontal and vertical motion is measured. These observations are then used to create particle size distributions (PSDs) and fall speed estimates at a one-minute resolution [16,19]. The images are subsampled by a factor of 25 to prevent repeatedly counting the same particle when calculating the PSD. The PIP camera images have an extremely low noise level, which enables a particle detection threshold of 0.3 mm equivalent area diameter. The PIP software makes computations for all diameter bins that contain at least one particle. Recent advancements in computer processing capabilities and expanded software since the development of the SVI enable the PIP to produce higher-order-derived data products, such as density and precipitation estimates and enable additional capabilities such as phase-separated rain and snow rate estimates [8,14,19,31]. These high-order data products are further described in Section 2.2.1. PIP-Derived Data Products. This work primarily focuses on the PIP measurements (e.g., PSDs, fall speeds, particle dimensions) and high-order-derived products.

2.1.2. Snow Field Observations

The NWS MQT office maintains an adjacent open field that is used for snow accumulation measurements. The seasonal snowpack depth is measured using a flat parcel of land 36 m² with four 150 cm snow stakes evenly distributed. New accumulations of snow are obtained five times per day by NWS MQT meteorologists at 0000, 0459 (one minute prior to local midnight), 0600, 1200, and 1800 UTC during ongoing precipitation events. A plastic, white snow board is placed adjacent to the open field and cleared off between each measurement [32]. The liquid water equivalent (LWE) precipitation is also collected and determined using a 20.32 cm diameter standard rain gauge. Snow from the gauge is melted to obtain the LWE accumulation and calculate the snow-to-liquid ratios (SLRs) for the 6-hourly periods (combining accumulation from the 0459 and 0600 UTC measurements). Observations from the snow field will be compared to the PIP-derived products and the retrievals in this work.

2.1.3. Accumulation Gauge

An OTT Pluvio 200 weighing precipitation gauge is also adjacent to the snow field. The Pluvio has an opening of a 16 cm diameter that is surrounded by a heated ring to prevent icing and snow build up. The Pluvio gauge is encircled by an Alter-type windshield of diameter 106 cm containing 24 lamellas to mitigate the impacts of wind on particle collection, thus improving catching efficiency of the gauge [31]. Accumulation observations are output at one-minute intervals with a resolution of 0.1 mm/min.

2.1.4. MicroRain Radar

The METEK Micro Rain Radar 2 (MRR) is a vertically profiling, 24-GHz frequency-modulated, continuous-wave Doppler radar [33]. The MRR has been used in field campaigns for observations of snow (i.e., [10,19,34,35,36]). The MRR at the NWS MQT site is configured to profile up to 3000 m above ground level (AGL) with a vertical resolution of 100 m and 1 min temporal resolution. The MRR dish contains a heater to keep the antenna free of snow and subsequent attenuation. However when snowfall is particularly dense, and/or the ambient temperature is near 0 °C, and/or the snow rate is very high, the melting from the heater can result in residual layers of water that attenuate the signal and impede the data. It is important to be aware of the impact of the heater during warm, heavy, and extremely dense snowfall events.

2.1.5. Surface Meteorological Observations

The surface meteorological data are obtained from the automated Davis Vantage Pro2 weather station, which is mounted on a tower (10 m AGL) adjacent to the NWS MQT office and snow field. Observations of surface temperature, dew point, and wind speed are recorded at one-minute intervals (instrument details can be found in Pettersen et al. [19]). The surface meteorological data are used to help choose the events used in this work, as low wind speed snow events will provide the best conditions for comparisons.

2.2. Methods for Estimating Snow Properties

Three methods of estimating bulk snow properties will be compared in this work. All three methods use the PSD and fall speed information from the PIP to produce an LWE snow rate (R):

$$R=\frac{1}{{\rho }_{liq}}\int m\left(D_e\right) \mathrm{V}\left(D_e\right) N\left(D_e\right) d{D}_e$$

(1)

where m(D_e) is the mass, V(D_e) is the fall speed, N(D_e) is the size distribution, and ${\rho }_{liq}$ is the density of water (note: all parameters discussed within this section are summarized in

Table A1. This table described the parameters for the PIP-derived data products in Section 2.2.1.

Parameter Symbol	Parameter Description	Units
De	area-equivalent diameter	mm
DVe	volume-equivalent diameter	mm
Dm	mass-weighted mean diameter	mm
${\rho }_e$	particle distribution equivalent density	g cm−3
$\overline{{\rho }_e}$	bulk equivalent density	g cm−3
${\rho }_{size}$	size-averaged density	g cm−3
${\rho }_{liq}$	water density	g cm−3
${\rho }_{min}$	minimum density boundary condition	g cm−3
R	liquid water equivalent precipitation rate	mm h−1
m(De)	mass distribution	g
V(De)	velocity distribution	m s−1
Vmax(De)	maximum fall speed for rain drop of De	m s−1
Vref	minimum fall speed boundary condition	m s−1
N(De)	size distribution	m−3 mm−1

in Appendix A). As part of the higher-order products created by the PIP processing software, the PIP method parameterizes particle density, herein referred to as equivalent density, to obtain the mass (see Section 2.2.1). The other two methods, von Lerber et al. [17] and Wood et al. [30], outlined in Section 2.2.2 and Section 2.2.3, respectively, retrieve mass directly. The von Lerber and Wood retrievals determine the parameters α and β of a power law, which describes particle mass as a function of the area-equivalent diameter (D_e) of the particle:

$$m\left(D_e\right)=\alpha D_e^{\beta }$$

(2)

Given the retrieved parameters and the PSD observed by the PIP, the results are then used to calculate mass-weighted mean diameter (D_m) and bulk equivalent density, the equivalent density averaged over the PSD. With information about particle fall speed as a function of size (observed by the PIP and handled differently by the two retrievals as described below), snowfall rates and corresponding accumulations can also be calculated.

2.2.1. PIP-Derived Data Products

As mentioned previously, the PIP provides higher-order data products estimating density and LWE precipitation rates. These products are derived using PIP observations of particle sizes and fall speeds in conjunction with an empirically determined parameterization of an equivalent density. Equivalent density, hereafter ${\rho }_e$, is the ratio of the mass of a nonrain particle to its apparent volume as determined from its area-equivalent diameter (D_e the diameter of a circle with the same area as that projected by a particle in the observed plane of view). The concept is similar to the snow particle bulk densities derived by Brandes et al. [37] and Huang et al. [38] from 2D Video Disdrometer (2DVD) imagery. The ${\rho }_e$ formulation is a data-guided parameterization that is determined by the D_e and size-resolved fall speed V(D_e) values observed by the PIP and is constrained by empirically derived physical limitations of ${\rho }_e$ and V(D_e). The ${\rho }_e$ parameterization includes assumptions about how density varies with particle size (D_e), and a fall speed-based interpolation between limiting values for density. The formulation can be written as:

$${\rho }_e(D_e, V\left(D_e\right))={\rho }_{liq}-\frac{{\rho }_{liq}-{\rho }_{min}\left({\mathrm{D}}_e\right)}{{\mathrm{V}}_{max}\left({\mathrm{D}}_e\right)-V_{ref}} \left({\mathrm{V}}_{max}\left(D_e\right)-V\left(D_e\right)\right).$$

(3)

The boundary condition of minimum density, ${\rho }_{min}\left(D_e\right)$, is derived based on previous work suggesting power law relationships for particle density. For example, Huang et al. [38] analyzed fall speeds of individual particles obtained from a 2DVD in Finland during light snowfall. They found that averaged size-resolved density (${\rho }_{size}$) as a function of volume-equivalent diameter (D_Ve; as opposed to the PIP use of area-equivalent diameter) was well represented by the following power law:

$${\rho }_{size}=0.15D_{Ve}^{-0.86}.$$

(4)

The PIP ${\rho }_e$ empirical parameterization boundary condition employs a similar power law model, as this approach is consistent with the use of power laws for particle mass and particle volume. The ${\rho }_e$ parameterization therefore uses the relation:

$${\rho }_{min}\left({\mathrm{D}}_e\right)=c{D}_e^x$$

(5)

where D_e is observed by the PIP, the gain (c) is based on empirical observations from a snow event (described below), and the exponent (x) is –0.86 following Huang et al. [38]. The boundary condition of the maximum density is simply the density of water, ${\rho }_{liq}$ (1 g cm⁻³).

The interpolation of density with fall speed of the ${\rho }_e$ parameterization is derived from physically based models for particle terminal velocity (e.g., Böhm [39]), which show, given a particular particle volume and shape, the particle terminal velocity increases as mass increases. This is also consistent with the particle density increasing, thus the PIP adopts a similar approach with which to extend the ${\rho }_e$ parameterization to higher-density (e.g., melting) precipitation events. If we assume for each D_e that the observed fall speed is the terminal velocity, as the density increases the terminal velocity increases up to the limiting case of the terminal velocities of raindrops. When precipitation falls as rain, the fall speeds should approximate those determined by Atlas and Ulbrich [40] terminal velocity of rain as a function of D_e. The Atlas and Ulbrich model is valid for raindrops up to 5 mm D_e. Therefore, the maximum fall speed as a function of particle size, ${\mathrm{V}}_{max}\left(D_e\right)$, is taken to be that of a raindrop of the same size and is estimated from expressions for raindrop fall speed by Atlas and Ulbrich [40]. The minimum boundary condition for the fall speed-based interpolation, V_ref, is empirically derived from PIP snow observations (detailed below).

The constants c, x, and V_ref were determined empirically from high-quality reference cases of snowfall event accumulations with data obtained from the PIP and snow field observations from the NWS MQT during the 2014–2015 winter. These events had SLRs ranging from low (10:1) to high (50:1) values, corresponding roughly to events with high and low snow particle densities, respectively. The snow event chosen to constrain the low-density end of the ${\rho }_e$ parameterization was a high SLR event with a mean ratio of 50:1, which had uniform snowfall and occurred under ideal meteorological conditions for accumulation accuracy (surface temperatures <−5 °C and wind speeds <5 m s⁻¹) to obtain estimates of ${\rho }_{min}\left(D_e\right)$. The event accumulation was determined using the following relationship:

$${\sum }_{event}{\left(R\Delta t\right)}_i={\sum }_{event} {[\frac{\pi }{6}\int c D_e^{-0.86} D_e^3 V\left(D_e\right) N\left(D_e\right)d{D}_e \Delta \mathrm{t}]}_i,$$

(6)

where ${\sum }_{event}{\left(R\Delta t\right)}_i$ is the total LWE accumulation summed over the duration of the event as determined by the NWS MQT snow field observations and the summation on the right is over all the discrete samples ($i$) observed by the PIP during the event. The right side of the equation is represented by variables obtained by the colocated PIP observations (D_e, the size-resolved velocity distribution $V\left(D_e\right)$, and the number concentration, $N\left(D_e\right)$). The gain (c) was then adjusted so that the resulting accumulation from the 50:1 SLR event matched that of the NWS MQT snow field (represented as ${\sum }_{event}{\left(R\Delta t\right)}_i$ in Equation (6)). The PIP observations of particle fall speeds during this high SLR (50:1) event also informed the minimum boundary condition for fall speed in the ${\rho }_e$ parameterization, resulting in the constant value for V_ref. The resulting ${\rho }_{e }$ parameterization was then tested and found to reproduce accumulations for both the low (10:1) and high (50:1) SLR NWS MQT observations. Additionally, the analyses in this work will further evaluate the performance of the ${\rho }_e$ parameterization more quantitatively.

The ${\rho }_e$ parameterization therefore allows for discrete, time-resolved precipitation rates for snowfall events, which are obtained from the PIP observations as:

$$R\left(t\right)=\frac{\pi }{6} \int \rho \left(D_e, V\left(D_e\right),c,x,V_{ref}\right) D_e^3 V\left(D_e\right) N\left(D_e\right) d{D}_e.$$

(7)

We have omitted the time dependence on the right side of the equation for clarity, and show explicitly the dependence of the ${\rho }_e$ parameterization on the described constants, as well as on observations of $D_e$ and $V\left(D_e\right)$. The three-dimensional shape of the surface of the ${\rho }_e$ parameterization as it relates to both D_e and fall speed is illustrated in Figure 1.

The PIP method uses the ${\rho }_e$ parameterization and observations of D_e and fall speeds ($V\left(D_e\right)$) to create a data product of ${\rho }_e$ (in g cm⁻³) resolved by particle size and time. Additionally, a one-minute bulk equivalent density, hereafter $\overline{{\rho }_e}$, is produced, which is the volume-weighted average of the ${\rho }_e$ distribution as shown in Equation (8):

$$\overline{{\rho }_e}=\frac{\frac{\pi }{6}\int {\rho }_e\left(D_e,\mathrm{V}\left({\mathrm{D}}_e\right)\right)D_e^{3}N\left(D_e\right)d{D}_e}{\frac{\pi }{6}\int D_e^{3}N\left(D_e\right)d{D}_e}$$

(8)

A single $\overline{{\rho }_e}$ value represents the equivalent density of all the snow that falls in the one-minute period. The value of $\overline{{\rho }_e}$ is approximately related to SLR. The ratio of $\overline{{\rho }_e}$ to the density of water (${\rho }_{liq}$) is unitless, and since the value of ${\rho }_{liq}$ is 1 g cm⁻³, the inverse of the $\overline{{\rho }_e}$ represents an SLR. Note that the equivalence of the PIP $\overline{{\rho }_e}$-derived SLR to the NWS-measured SLR depends on the PIP acquired volume accumulation being consistent with the resulting volume of the snow as it lies on the ground [41]. It is important to note that the ${\rho }_e$ parameterization was developed and tested on only a few snow events with SLRs ranging from 10:1 to 50:1 (the limiting case). The presented work examines a larger sample size of snow events and includes a broader range of SLRs and will assess the ${\rho }_e$ parameterization for snow conditions outside of the initial training events. The snowfall events selected for comparison in this work have a range of NWS SLRs from 7:1 (wet, dense snow) to 70:1 (ultralow density) and therefore span $\overline{{\rho }_e}$ values of 0.14 to 0.014 g cm⁻³, respectively. These values are typical for precipitation events at MQT that are classified as snow, but do not approach the expected ${\rho }_e$ values for moderately melted snow, as would be found in mixed-phase precipitation, and for rain, which is an area for future study. In this work, we will compare the PIP-derived data products of the bulk density ($\overline{{\rho }_e}$) and LWE precipitation rate against results from two previously established snow mass retrievals: von Lerber et al. [17] and Wood et al. [30] described in the following sections, as well as compared to observations from the NWS MQT snow field 6-hourly accumulation measurement.

2.2.2. von Lerber Mass Retrieval

The mass retrieval process, described in more detail in von Lerber et al. [17], is based on the hydrodynamic theory [39], in which, the gravitational force related to the particle mass is equalized with the atmospheric drag influencing the falling particle. In this method, the mass is retrieved for a single particle utilizing the observations of fall velocity, dimensions, and environmental variables, such as the temperature, relative humidity (obtained from the ERA 5 reanalysis), and pressure. In the literature, there exist several parameterizations to describe particle flow characteristics (i.e., the relation between the Best and Reynolds numbers), and in the method used here four of these parameterizations are computed: Böhm [39], Mitchell and Heymsfield [42], Khvorostyanov and Curry [43], Heymsfield and Westbrook [44]. The parameterization that best represents the estimated snow accumulation is selected for each individual snow event. The choice is based on scaling performed by comparing the estimated total LWE accumulation of the mass retrieval process with a total observed accumulation for each snowfall event. Whenever a significant difference between the estimated and observed accumulation is seen, a simple correction factor, ranging from 1.0 to 1.42 is applied to the particle maximum diameter (diameter of the circumscribing sphere) when calculating the Reynolds number. Typically, this factor has a value of 1.0 and no correction is performed. The idea behind this correction factor is to balance the error in the irregular-shaped particle dimensions introduced by the observation from a single projection plane perpendicular to the fall direction [7,45]. In addition, this correction factor is also seen implicitly to compensate truncation of the observed PSD [17]. Thus, the uncertainties of this mass retrieval process are mostly related to observations of the dimensions (as they are performed from the side and not from the fall direction), PSD truncation in observations, selection of the used parameterization, and impacts from gauge undercatchment used for scaling. In this study, the scaling is performed against the Pluvio gauge adjusted by the NWS MQT snow field observations.

Prior to the scaling, between the estimated and observed LWE accumulation, the PIP obtained one-minute PSDs are filtered. First, by applying a threshold in a total number concentration of 20 m⁻³; second, by removing the isolated observations of large particles (excluding outliers if four smaller size bins are not containing any particles, i.e., Tiira et al. [15]); and third, by taking a 5 min mean. The coefficients of $m\left({\mathrm{D}}_e\right)$ in power-law are derived through a linear regression fit to the 5 min retrieved masses in log-scale, requiring at least 30 observed particles, and outliers are filtered utilizing Gaussian kernel density estimation (KDE; [46]) to find the most probable mass for each diameter bin, and considering only mass observations within half-width at half maximum from the bin peak KDE value [15]. The $\overline{{\rho }_e}$ is then determined every 5 min as a function of ${\mathrm{D}}_e$ with the coefficients of $m\left({\mathrm{D}}_e\right)$ and the averaged observed PSD.

2.2.3. Wood Mass Retrieval

To leverage the information provided by multiple sources of snowfall observation during intensive precipitation ground validation experiments, Wood et al. [13,30] devised an optimal estimation method to retrieve the parameters described above for snow particle mass (Equation (2)), as well as the parameters of a power law describing horizontally projected area and a shape parameter relating the observed particle maximum dimension to its likely true dimension. The retrieval uses measurements of snow PSD, radar reflectivity, particle fall speeds, and precipitation rate along with an a priori probability density function for the parameters to be retrieved. This combination of observations was found to effectively constrain the mass parameters and provide less robust constraints on the area and shape parameters [30]. Given the mass and projected areas, fall speeds are forward-modeled using hydrodynamic balance [39] as implemented by Mitchell and Heymsfield [42]. The retrieved shape parameter compensates for differences between the PIP-observed particle size and the actual particle maximum dimensions [7]. As diagnostics, the retrieval fits the observed PSDs to a negative exponential function using orthogonal distance regression [47] to obtain PSD slopes and intercepts.

The PIP provides both PSDs and size-resolved fall speeds, reflectivities are from the MRR, and precipitation rates are from the combined Pluvio and snow stake field observations as described in Section 2.1. As in Wood et al. [13,30], the size-resolved falls speeds are composed into mean fall speeds over three size ranges, which for this work were set to D_e values of 2.75 ± 0.25 mm, 1.5 ± 0.125 mm, and 0.75 ± 0.125 mm. These ranges are intended to provide information about how fall speed varies with size without overemphasizing the importance of the fall speed information in the retrieval. If no fall speeds are observed in one or more of the size ranges, the retrieval is not performed for that time step. For this work, 1261 retrievals out of 4632 samples were not performed, leaving 3371 samples for retrieval. This was most frequently due to observing no particles in the largest size range. Radar reflectivities are taken from the fourth range bin above the surface, approximately 120 m above the surface, to avoid contamination by ground clutter. As in Wood et al. [13,30], observed, not fitted, PSDs are used for the retrieval.

Uncertainties in observed PSD, observed and forward-modeled fall speed, and observed and forward-modeled precipitation rate follow Wood et al. [7,13,30]. For reflectivity at the MRR frequency, the combined model-measurement uncertainties involve not only measurement errors, but also model errors due primarily to unknown particle shape. These uncertainties were modeled assuming a constant 1 dB measurement uncertainty and a reflectivity-varying model uncertainty based on an analysis of a higher-frequency (W-band) radar [48], which would exhibit somewhat more shape sensitivity compared to the MRR. During light snowfall (high SLR) conditions, the variances for the reflectivity model-measurement were inflated by a factor of 2 to account for the vertical separation between the observations. During heavy snow (low SLR), a factor of 40 was used so that the measurements from the Pluvio were emphasized by the retrieval.

The nature of the retrieval requires reasonable physical consistency among the observations. If not, the retrieval may fail to converge or converge to unphysical results. For these observations, the principal concern is the relatively weak sensitivity of the Pluvio to lighter snowfall compared to the PIP and snow field observations. The most common manifestation of this inconsistency involved samples for which the PIP detected nontrivial concentrations of snow particles and the radar reflectivity was consistent with precipitation, but the Pluvio detected no accumulation during the sample. The MRR reflectivity was often substantial for these samples, with a mean value of 8.2 dBZ (figure not shown). For such samples, the retrieval typically gave results that produced a small snowfall accumulation well below the 0.1 mm detection limit of the Pluvio. Though some of these samples may represent valid retrievals during very light snowfall, a conservative filter was applied to the output that removed any samples in which the Pluvio reported no accumulation. The Pluvio failed to detect precipitation in 2062 samples out of the 3371 samples for retrieval, with most of the failures (1561) occurring during the high SLR events.

2.3. Snow Event Selection

In selecting the snow events to use in this work, we purposely set thresholds to limit uncertainties originating from environmental conditions. Accurate comparisons of bulk snow characteristics, such as LWE precipitation rate, mass, and density, can be impeded by environmental factors, such as moderate to high wind speeds and surface temperatures nearing the melting point. We therefore set limits of surface temperatures <0 °C (slightly more conservative than the <2 °C limit that were used in previous studies of snowfall at the site [19]) and wind speeds <5 m s⁻¹ for the event duration. Additionally, we selected snow events that had continuous snowfall for at least six hours, so that we could compare the outputs to at least one snow field accumulation measurement. Finally, we partitioned the events into low versus high SLR snow categories by using the 6-hourly NWS MQT SLR observations. By using the NWS observed SLR values to categorize the snow events, we partition only on differences in the density of snow on the ground and not by information about the particle densities, in which number density, size, and riming will all have impacts [49]. The low SLR snow events are those where the SLR is 15:1 or less, and the high SLR snow events are where the SLR is greater than 15:1. Using these criteria, we acquired 11 low and 8 high SLR snow events from the NWS MQT site between January 2017 and April 2019.

3. Results and Discussion

3.1. Low and High SLR Snow Event Characteristics

The low SLR snow events are shown in

Table 1. This is a table of the 11 low SLR snow events used in this work. Listed are the start and end times used for the analyses, as well as mean values of snow-to-liquid ratio (SLR), snow rate (here abbreviated as SR), wind speed (here abbreviated as WS), and temperature for the duration of each event.

Start (Date/Time)	End (Date/Time)	Mean SLR NWS	Mean SR PIP (mm h−1)	Mean WS (m s−1)	Mean Temp (°C)
11 November 2017 0600 UTC	11 November 2017 1200 UTC	12.14	0.899	3.33	−5.25
11 December 2017 1300 UTC	11 December 2017 1900 UTC	15	0.095	3.28	−8.95
13 December 2017 1000 UTC	13 December 2017 1800 UTC	13	0.278	1.55	−11.49
4 January 2018 0000 UTC	6 January 2018 1600 UTC	15.03	0.137	3.42	−17.16
22 January 2018 1800 UTC	23 January 2018 0300 UTC	7.5	3.751	3.96	−3.53
3 February 2018 1800 UTC	4 February 2018 0400 UTC	12	0.316	1.83	−12.88
12 April 2018 0800 UTC	12 April 2018 1600 UTC	10.45	1.672	0.88	−0.46
15 April 2018 1000 UTC	15 April 2018 2300 UTC	7.34	3.134	6.00	−6.41
2 December 2018 0600 UTC	3 December 2018 0300 UTC	9.11	2.455	4.73	−1.74
7 January 2019 0700 UTC	7 January 2019 2300 UTC	12.65	1.035	3.47	−2.22
20 January 2019 0600 UTC	21 January 2019 0600 UTC	14.89	0.274	2.87	−16.89

. These snow events have mean SLR values ranging from approximately 7.5:1 (extremely wet, dense snow) to 15:1 (commonly specified relationship for synoptically forced snowfall [50,51]). The mean PIP-derived snow rates vary widely from 0.095 to 3.75 mm h⁻¹. The observed event mean surface temperatures also have a large spread, ranging from −17 to −0.5 °C. There is a tendency for the warmer low SLR events (those with mean temperatures closer to 0 °C) to have higher mean snow rates. Generally speaking, the majority of these low SLR snow events are produced by synoptically forced, deep precipitation (figures not shown). However, two of the low SLR events would be categorized as shallow LeS: 4 January through 6 January 2018 and 20 January through 21 January 2019. Both of these low SLR LeS events have extremely cold temperatures (event means of −17 °C), which likely indicates that boundary layer temperatures were colder than the dendritic growth zone (DGZ) and would not be an environment conducive to large particle growth [52]. This is consistent with the characteristics of a secondary mode of LeS snow identified at MQT in previous work [19,20].

Characteristics of the high SLR snow events are shown in

Table 2. This is a table of the eight high SLR snow events used in this work. Organization follows that of Table 1.

Start (Date/Time)	End (Date/Time)	Mean SLR NWS	Mean SR PIP (mm h−1)	Mean WS (m s−1)	Mean Temp (°C)
14 December 2017 0000 UTC	14 December 2017 0800 UTC	32.75	0.259	2.59	−10.21
15 January 2018 2200 UTC	16 January 2018 2000 UTC	36.39	0.351	1.00	−9.76
28 January 2018 0800 UTC	28 January 2018 2300 UTC	70	0.054	2.48	−8.86
6 March 2018 1800 UTC	7 March 2018 2300 UTC	27.58	0.511	3.24	−7.26
13 March 2018 0400 UTC	7 March 2018 2300 UTC	30.33	0.132	3.94	−7.78
17 November 2018 0000 UTC	17 November 2018 2100 UTC	29.62	0.148	3.14	−6.04
25 November 2018 1200 UTC	26 November 2018 2000 UTC	30	0.079	3.49	−5.43
10 January 2019 0000 UTC	10 January 2019 2300 UTC	59	0.053	3.39	−10.50

. The high SLR snow events have mean SLRs from 28:1 (lower density snowfall) to greater than 50:1 (very low density), which are consistent with SLRs seen in previous studies of LeS [20,53,54]. Originally, we identified two additional ultralow density snow cases (SLR > 80:1); however, the Pluvio did not record any data during these events, which makes it difficult for the von Lerber and Wood methods to accurately retrieve mass. The mean PIP-derived snow rates are much lower than the low SLR events and a narrower range of values, from 0.079 to 0.511 mm h⁻¹. The event mean surface temperatures ranged from −10 to −5 °C, again a much smaller distribution than that of the low SLR snow events. All of the high SLR events were a product of surface and air interactions, namely LeS, orographic forcing, or both, as supported by the precipitation radar profiles (figures not shown).

3.2. Snow Microphysical Characteristics

The PIP produces end-user products such as PSDs and size-resolved particle fall speeds, which illustrate microphysical characteristics of the snow. Figure 2 shows the PIP-derived PSDs for the low and high SLR events presented as two-dimensional histograms (Figure 2a,c, respectively). The low SLR snow events have large particle concentrations (>10³ m⁻³ mm⁻¹) at small D_e (<1 mm), while the high SLR snow events demonstrate a broader range of PSDs with an order of magnitude fewer small particles (most values <10³ m⁻³ mm⁻¹). These PSD results are consistent with prior studies of snowfall regimes [19]. The composite fall speed distributions also indicate differences between the low and high SLR events (Figure 2b,d, respectively). The low SLR snow composite implies generally higher fall speeds as a function of D_e compared to the high SLR snow composite, which has a flatter relationship and smaller particle fall speeds (<1 m s⁻¹). The mean fall speeds for the low SLR snow events are 30–45% larger than those for the high SLR snow events as a function of D_e for those D_e values <5 mm.

We also can examine the microphysical characteristics of the low and high SLR cases through differences between the shape of the size distributions for each category. Previous studies have applied the following inverse exponential relation to measured snow particle concentrations:

$$N\left(D_e\right)=N_0 \exp \left(-\mathsf{\Lambda } D_e\right)$$

(9)

where $N\left(D_e\right)$ is the particle concentration per particle size, D_e. N₀ and Λ are the PSD intercept and the slope parameter, respectively [19,55,56,57,58]. Both the von Lerber and Wood retrieval methods applied the relation in Equation (9) to 5 min intervals of the PIP data and each obtained N₀ and Λ for the low and high SLR events. There are differences in the resulting N₀ and Λ derived by the von Lerber and Wood methods due to differences in fitting methods and preprocessing of the PIP PSD data (see Section 2.2.2 and Section 2.2.3 for details).

Figure 3 illustrates composites of the derived von Lerber and Wood N₀ and Λ for the low SLR events (Figure 3a,b, respectively) and the high SLR events (Figure 3c,d, respectively). In general, the low SLR composites from the two methods show similar ranges, with N₀ ranging from 10² to 10⁵ m⁻³ mm⁻¹ and Λ values of 0.5 to 6.5 mm⁻¹ in a seemingly linear structure (see Figure 3a,b). The von Lerber shows fewer small values of both N₀ and Λ, which could be due to preprocessing that applies the threshold for total concentration and removes single occurring particle counts prior to the fitting routine (see Section 2.2.2). In both the von Lerber and Wood outputs, the high SLR events N₀ and Λ values are much smaller than that for the low SLR events, with N₀ ranging from 10 to 10⁴ m⁻³ mm⁻¹ and Λ ranging from 0.4 to 1.6 mm⁻¹, and no apparent relationship between the intercept and slope parameter. These results clearly show there are different shapes to the PSDs as a function of events that are low versus high SLR, which align with the results illustrated in Figure 2a,c. The median von Lerber values of both N₀ and Λ are 66% smaller for the high versus the low SLR snow events, while the median Wood value of N₀ is 58% smaller and Λ is 49% smaller for the high versus low SLR snow events. The N₀ and Λ results from the low SLR composite imply that these snowfall events have narrow PSDs with large amounts of small particles, similar to what has been observed for previous synoptically forced snow [19,55,58,59,60]. While the high SLR composite indicates the opposite: that these snow events tend to have broad PSDs with an order of magnitude fewer small particles and a tendency towards containing more large particles. The N₀ and Λ are similar to values seen in previous studies of LeS [19,56,61].

Additionally, we examined and compared the retrieved values of mass-weighted mean diameter for each SLR snow regime and the von Lerber and Wood methodologies. Mass-weighted mean diameter (D_m) is:

$$D_m=\frac{{\int }_{D_{min}}^{D_{max}}m\left(D_e\right)N\left(D_e\right)D_{e}dD}{{\int }_{D_{min}}^{D_{max}}m\left(D_e\right)N\left(D_e\right)dD}$$

(10)

There is good agreement (correlation coefficients (cc) of 0.945 and 0.883) between the von Lerber and Wood-derived values for D_m (see Figure 4). The high SLR event composite of D_m indicates a larger spread between the von Lerber and Wood values versus the low SLR, with the Wood D_m values generally higher for both categories (biases of +13.7% and +10.1%, respectively). In the low SLR event composite, both methods indicate that the majority of D_m values are smaller than 2 mm. However, in the case of the high SLR events, the range of values of D_m is much larger, extending from 1 to 4 mm, with a concentration towards larger values of D_m mostly between 1.5 and 3 mm. The median von Lerber value of D_m is 51% larger and the median Wood value of D_m is 56% larger for the high versus the low SLR snow events. The GPM Dual Precipitation Radar (DPR) algorithm has prescribed relationships for precipitation and D_m for stratiform and convective rainfall [21,62], and it is important to also examine these relationships for snowfall. Figure 4 illustrates that there are clear D_m differences for the low versus high SLR snowfall events. The SLR event types occur under different environmental conditions that should be therefore an additional factor that is considered when examining precipitation and D_m relations.

In addition, we investigated the output of the mass parameters, α and β, described in Equation (2), from the von Lerber retrieval (Figure 5a,c) and the Wood retrieval (Figure 5b,d). Note that the Wood retrieval outputs are fewer due to Pluvio limitations and filtering detailed in Section 2.2.3. (especially true for the high SLR cases). The low SLR events have scaling factor values, α, slightly larger than 10⁻⁴ and exponent, β, between 2 and 2.5 as produced by the von Lerber retrieval, while the Wood retrieval has a lower range of α (between 10^−4.25 and 10⁻⁴) and similar β (2–2.5). The von Lerber high SLR event composite shows much lower values for α (between 10^−4.75 and 10⁻⁴) and β (between 1.75 and 2.25). Though there are few points, the Wood retrieval also indicates lower α values, whereas the β range is similar to the low SLR output. Similar to the N₀ and Λ comparisons, we see that the low and high SLR snow event composite relations for α and β exhibit a clear separation, as illustrated in the von Lerber retrievals (Figure 5a,c). The values of both α and β are larger for the snow particles observed by the PIP during the low SLR events, which indicates that the mass increases faster as the area-equivalent diameter of the particle increases compared to the high SLR particles. The median of the von Lerber α values for the low SLR snow events is twice that of the high SLR, while the medians of the von Lerber β values are approximately the same for the low and high SLR (difference of 6%). The differences seen between the low versus high SLR events imply that snowfall retrieval algorithms, which rely on a priori assumptions of α and β either explicitly [48] or implicitly (e.g., Kulie and Bennartz [63], Liu [4], and Table 1, Braham et al. [64]), might benefit by adjusting to variations in snow type.

Finally, we compare the ice water content (IWC) retrieved by both the von Lerber and the Wood methods, for the low and high SLR snow events (see Figure 6a,b, respectively). Figure 6a evaluates the respective retrieved IWC values for the low SLR events. Though the IWC retrieved values are correlated (cc of 0.864) for the low SLR snow events, the Wood retrieval calculates much smaller IWC values versus the von Lerber retrieval, with a bias of –50%, and spread of 0.08 g m⁻³ as determined by the Root Mean Square Error (RMSE). In general, low SLR snow events produce more IWC values below 0.6 g m⁻³ for the von Lerber and 0.4 g m⁻³ for the Wood retrievals, respectively. The high SLR event composite has far fewer values to compare and there is a wide spread in the retrieval comparison and reduced correlation (cc of 0.539). There is an indication that the Wood output retrieves larger values of IWC matched to that of the von Lerber output; however, with so few points we can only suggest this is a possibility. There is evidence that the low SLR events produce more integrated mass as the median von-Lerber-retrieved value of IWC is more than an order of magnitude larger for the low SLR versus the high SLR snow events.

In general, the low SLR snow events have an order of magnitude more small particles and 30 to 45% higher fall speeds than the high SLR snow events. Both the von Lerber and Wood-derived values of N₀ and Λ are smaller (66% and ~50%, respectively) for the high versus low SLR snow events, which is consistent with findings from Pettersen et al. [19]. The D_m values for the high SLR snow events are larger for von Lerber (51%) and Wood (56%) versus the low SLR events. The median von-Lerber-retrieved value of α was two times larger for the low SLR versus the high SLR snow events, while the β values were approximately the same. The von-Lerber-retrieved IWC was more than an order of magnitude larger for the low versus the high SLR snow events. In terms of comparisons between the retrievals, the Wood D_m values were between 10% and 14% larger than the von Lerber D_m values, regardless of snow event type, and the D_m values for the retrievals were highly correlated (cc values of 0.88 and 0.95). The Wood-retrieved values of IWC were biased 50% lower than the von Lerber IWC values for the low SLR snow events; however, the retrieved values were highly correlated between the retrieval methods (cc of 0.86).

3.3. Bulk Snowfall Characteristics

We have shown that using snow observations from the PIP the von Lerber and Wood-retrieved microphysical characteristics compared well and showed consistent differences between the low and high SLR events. Additionally, the values of α, β, and mass were physically reasonable and again illustrated distinct relationships for the low and high SLR snow events. Here, we examine and intercompare the bulk characteristics of the low and high SLR snow events from the PIP higher-order-derived data products and the von Lerber and Wood retrievals. The intercomparison focuses on the produced values of snow rate and $\overline{{\rho }_e}$.

Figure 7 highlights comparisons of the PIP-produced LWE precipitation rates (on the y-axes) with those from the von Lerber (Figure 7a,c) and Wood (Figure 7b,d) methods for the low (Figure 7a,b) and high (Figure 7c,d) SLR snow events. For the low SLR event composites, we see LWE precipitation rates up to 6 mm h⁻¹ with a very high correlation between the PIP-derived snow rates and the output from von Lerber (cc of 0.976), while less correlated with the snow rates retrieved from Wood (cc of 0.786). In general, the LWE precipitation rates are less than 2 mm h⁻¹, but there are indications of repeated instances of higher rate values. The PIP LWE precipitation rates tend to be in good agreement with those of von Lerber, a bias of +1.56%, but biased high with respect to Wood (+7.41%) for the low SLR snow events. For high SLR snow events, the PIP and von Lerber retrievals compare well (cc of 0.940), with a linear relationship as LWE precipitation rates increase. There is again good agreement between the produced snow rates, with the PIP is biased slightly high (+3.1%) compared to the von Lerber values. Many of the high SLR snow event samples provided physically inconsistent inputs for the Wood method and were so excluded (see Section 2.2.3), so it is difficult to compare. However, we can compare the differences between the low and high SLR snow events by examining the PIP and von Lerber produced values of snow rate. Though the RMSE between the PIP and von Lerber values of snow rate is 0.115 mm h⁻¹ for the low SLR and 0.242 mm h⁻¹ for the high SLR snow events, the median values for the event regimes agree very well, with a 4% difference for the low SLR and a 9% difference for the high SLR snow events. Unsurprisingly, the low SLR snow events have much larger LWE precipitation rates, with median rates more than 20 times higher compared to the high SLR snow events for both the PIP and von Lerber outputs.

We also assess the respective values of $\overline{{\rho }_e}$ (g cm⁻³) for the low and high SLR snow events (see Figure 8). As discussed in Section 2.2.1, the SLR output from the three methods is roughly the inverse of the ratio $\overline{{\rho }_e}$ to the density of water (1 g cm⁻³); however, it is important to note that this conversion does not account for compaction effects. Here, we compare the $\overline{{\rho }_e}$ values by inverting them to produce bulk SLR values for each time step. The low and high SLR events were defined using SLR observations from the NWS MQT snow field, and therefore we expect the comparisons in Figure 8 to have different ranges of PIP, von Lerber, and Wood produced SLRs as a function of event category, regardless of method. We see for the low SLR snow event composites (Figure 8a,b) that the vast majority of SLR values are less than 10:1 for the PIP, von Lerber, and Wood retrievals. The SLR output from the PIP and the von Lerber retrievals have a correlation of 0.756. There is a shift for SLR of >15 where von Lerber reports larger values compared to the PIP, which results in the PIP values biased 30% lower than output from von Lerber. However, the median values of the PIP and von-Lerber-derived SLR are within 6%. The low SLR composite for the PIP and Wood-derived SLR have a low correlation (cc of 0.382), and an inclination for the PIP SLR to be much lower than that of the Wood output (bias of –79.6%). The high SLR event composite between the PIP and von Lerber SLRs (Figure 8c) show a similar correlation to the low SLR composite (cc of 0.804); however, with more variability (RMSE of 4.843), and a low bias of almost 50%. Additionally, the median PIP-produced SLR value is 60% smaller than that of von Lerber for the high SLR snow events. There are too few values of the Wood method for the high SLR snow events to conclude similarities or differences with the PIP output (Figure 8d).

The PIP and von Lerber low (Figure 8a) and high (Figure 8c) SLR event composites illustrate that the majority of the prescribed low SLR snow event values are <10:1, whereas the high SLR snow event SLRs range from 8:1 to 30:1, with some of the von Lerber values extending to 60:1 (note different x and y limits). Additionally, the median PIP value for the high SLR snow events is twice that of the low SLR snow events, while the corresponding von Lerber median values are 3.6 times larger. The range of SLRs values for the high SLR snow events is much larger than that of the low SLR snow events.

Though much of the focus of this work is an intercomparison of the bulk snow products from the PIP higher-order products and the von Lerber and Wood retrievals, we additionally wanted to assess each method against an established means of acquiring accumulation to produce uncertainty estimates. Here, we examine the six-hourly accumulation and SLR values from NWS MQT snow field measurements for all the presented snow events against corresponding intervals of accumulation. These accumulation comparisons are presented in Figure 9 with the NWS MQT snow accumulation observations on the respective x-axes (left panels). There is a high correlation between the six-hourly accumulations from the respective methods and the snow field observations, with PIP cc of 0.971, von Lerber cc of 0.982, and Wood cc of 0.986. The PIP-produced and von-Lerber-retrieved accumulations show excellent agreement with the snow field observations, where each is biased slightly high at +2.48% and +7.22%, respectively, The Wood-retrieved accumulations also agree very well; however, are biased slightly lower than the snow field at −9.45%. The accumulations for all three methods RMSE range from 0.11 to 0.13 mm per 6 h interval. The PIP (Figure 9a) and von-Lerber-derived (Figure 9c) accumulation comparisons have 46 and 44 6 h interval matches, respectively, while the Wood accumulation only had 22 6 h accumulation comparison points due to the Pluvio filtering method (see Section 2.2.3).

Additionally, we compare the six-hourly observations of SLR from the NWS MQT snow field (x-axes in Figure 9b,d,f) and the median SLR produced by each method. The PIP, von Lerber, and Wood median SLR tended to be lower than that measured by the snow field with biases of −68.2%, −56.4%, and −41.7%, respectively; however, all are well correlated (cc of 0.831, 0.874, and 0.843, respectively). The high SLR snow event median SLRs (purple circles) tend to have much more spread, which is unsurprising given the large variability in LeS and orographic snow events. There is less agreement between the snow field and the Wood SLR values; however, there are also fewer points (22) to compare. In general, the Wood median SLRs for the low SLR events agree better with the snow field than the retrieved values for the high SLR events. Additionally, the Wood median SLRs for the same time periods as the PIP and von Lerber median SLRs correspond well for the low SLR event samples; however, they do not for the high SLR events (figure not shown).

The results imply fundamental differences in the bulk snow characteristics between the low and high SLR snow regimes. The low SLR snow events produced much larger LWE snow rates, as both the PIP and von Lerber indicated 20 times higher values compared to the high SLR snow events. The differences in snow rates seen between the low SLR (primarily synoptically forced snow events) and high SLR (LeS and orographically forced snow events) are consistent with differences observed in previous studies [19,20,37]. Unsurprisingly, the PIP and von Lerber determined median values of SLR were 2 and 3.6 times larger, respectively, for the high versus the low SLR snow events. The high SLR snow events tended to have a much larger spread of values of PIP and von Lerber determined SLR, which could be due to the large variability inherent in surface-based processes by different forcing mechanisms (LeS versus orographic) and boundary-layer conditions [19,20]. In terms of intercomparisons of snow rates, there was excellent agreement between the three methods for the low SLR snow events, with the PIP values biased +1.58% compared to the von Lerber retrieval and +7.41% compared to the Wood retrieval. For the high SLR snow events, there was again excellent agreement with the PIP values biased +3.1% compared to the von-Lerber-retrieved values. When compared to the NWS snow field measurements of six-hourly accumulation, the retrievals all showed consistent results: The PIP accumulations were only biased +2.48% higher with an RMSE of 0.13 mm per 6 h time interval; the von Lerber retrievals of accumulation had a high bias of +7.22% and RMSE of 0.11 mm; and the Wood retrievals were biased −9.45% lower than the snow field accumulation with an RMSE of 0.13 mm per 6 h period. It is important to note that the PIP-produced accumulations agreed best with the NWS snow field and, unlike the Wood and von Lerber retrievals, do not use any information from the snow field to calculate the event accumulations. Although the scatter shown in Figure 7 suggests the fractional uncertainties in high-time resolution snow rates have larger spread (especially when compared to the Wood retrievals), Figure 9 indicates that 6 h PIP estimates of accumulation or mean snowfall rate appear to have uncertainties smaller than 5%.

3.4. Detailed Example of a Snow Transition Event

Throughout this analysis, we have presented the output from the PIP, von Lerber, and Wood methods using a composited approach that highlights the differences between the low and high SLR snow events, while simultaneously showing that the PIP higher-order-derived products accurately represent the bulk snow conditions. We additionally want to emphasize the utility of these methods to examine the snow characteristics as a function of individual events. In this subsection, we present a snow event that transitions from low to high SLR snowfall under changing environmental conditions. Detailed examples of low and high SLR events can be found in Appendix B (Figure A1 and Figure A2, respectively). The event is presented with ancillary data from the profiling radar (MRR), the Pluvio, and snow field accumulations, the surface meteorological conditions, and the output from the von Lerber and Wood retrievals to give further context and support to the PIP-derived data products.

We examine an event that started as synoptically forced snow and transitioned to lake-enhanced snowfall [20,65] from 19 November (0400 UTC) to 20 November 2018 (0200 UTC), shown in Figure 10. During the beginning of the event (prior to 1300 UTC), the MRR profile indicates a deep precipitation column of light snowfall (reflectivities <10 dBZ; Figure 10a) and comparatively slow fall speeds (<1 m s⁻¹; Figure 10b). We also see that the PIP PSDs are very narrow during this time period with particles smaller than 5 mm D_e (Figure 10c). Additionally, the PIP ${\rho }_e$ distribution values are greater than 0.1 g cm³ (10:1 SLR) during the early part of the snow event (Figure 10d). At 1300 UTC, we see an abrupt change in the observations with increasing reflectivity and Doppler velocity values with decreasing precipitation height (~1.5 km AGL). At this same time, the PIP PSDs respond with an immediate broadening (particles ranging up to 15 mm D_e), while still maintaining a high concentration of small particles (>10⁴ particles smaller than 1 mm D_e). We also see that the ${\rho }_e$ distribution shifts to lower values, generally less than 0.1 g cm⁻³ (Figure 10d). The accumulation (Figure 10e) is gradual during the low SLR period of the snow event, but the snow rate increases directly following the transition at 1300 UTC, with ~9 mm accumulation in 6 h (in agreement with the von Lerber accumulation). Overall, the von Lerber total accumulation exactly matched the NWS MQT snow field, while the PIP recorded a relatively higher amount (120%) and the Wood retrieval much less (67%). Unsurprisingly, we see an instantaneous response from the values of SLR and D_m at 1300 UTC as well (Figure 10f,g, respectively), and both increase rapidly, with SLR values from 10:1 to greater than 20:1 and average D_m growing from ~0.5 to more than 2 mm. During the transition period, winds shift to slightly higher speeds (Figure 10h) and more northerly direction.

4. Conclusions

The primary goal of this work was to define and demonstrate the quality and efficacy of the PIP higher-order products. The results detailed in this work confirm the ability of the PIP-derived products to quantify properties of LWE snow rate, ${\rho }_e$, and SLR. Specifically, through comparisons with established mass retrieval methods from von Lerber et al. [17] and Wood et al. [30], we have shown that the PIP produces physically consistent bulk snow characteristics and produces LWE snow rates that agree with these methods within 1.5 to 12% depending on snow type and retrieval. Additionally, the PIP-derived LWE accumulations over 6 h periods match the NWS MQT snow field measurements for both the low and high SLR events, suggesting that the PIP LWE snow rates are accurate for the examined snowfall. The PIP accumulations were within 2.48% of those for the NWS snow field observations and, unlike the von Lerber and Wood retrievals which made use of the snow stake field observations directly or indirectly, were produced with no ancillary information of precipitation rate or accumulation. The von Lerber and Wood-derived six-hourly accumulations also showed consistent results with the NWS MQT snow field measurements, agreeing within 10%. Furthermore, the results from this work indicate that the PIP is much better at detecting small accumulations of snow compared to the unfenced Pluvio gauge. The large number of samples for which the Pluvio detected no accumulations and the PIP (and MRR) indicated significant snowfall, particularly prevalent for the high SLR events, underscores this point. These results imply that the PIP is less impacted by wind and therefore has less or no undercatch. This may be due to the advantage that the PIP is able to sample in a large and unimpeded volume, where the Pluvio has only a small opening for which to receive falling particles. The PIP is therefore well suited to providing time series of snowfall rates and SLR at a high-time resolution, as is exhibited in the detailed analysis of specific snow events (see Figure 10, Figure A1, and Figure A2).

We also found that there are distinct microphysical and bulk characteristics for the low and high SLR snow events, as illustrated by the three retrieval methods. The PIP-produced PSDs, fall speeds, and relationships of N₀ and Λ were very different for the low versus high SLR and consistent with values found in previous studies of synoptically forced and lake-effect snowfall. In general, the D_m values for the high SLR snow events are larger and more variable than those for the low SLR snow events. Similar to the N₀ and Λ comparisons, we see that the low and high SLR snow event composite relations for α and β exhibit a clear separation, with lower values of both for the high SLR snow events. There were too few IWC values for the high SLR snow events to compare the von Lerber and Wood methods, but low SLR IWC values are generally larger. In general, the PIP LWE snow rates compared well to those produced by the retrieval methods for the high and low (von Lerber only) density events. The PIP and von-Lerber-retrieved values of SLR were consistent for both low and high SLR snow events, with the von Lerber method producing consistently larger values. Finally, accumulation comparisons between all three methods and the six-hourly snow field were excellent for both the low and high SLR snow events.

Future work will focus on a few key areas. First, to develop a means to produce real-time uncertainty estimates for the PIP LWE snow rate and ${\rho }_e$ products. Second, to enhance the current PIP ${\rho }_e$ parameterization. As outlined in Section 2.2.1, the PIP ${\rho }_e$ parameterization is currently constrained and tested by empirical observations from only three events with a limited range of SLRs. Further analysis of extreme high SLR, higher density (<10:1), and mixed-phase snow events, as well as those containing a majority of graupel and rimed particles will help to improve the ${\rho }_e$ parameterization. Third, continue to expand our understanding of the meteorological and environmental factors on the retrieved snow micro- and macrophysical characteristics. This work indicates that there are specific microphysical relations and bulk characteristics that are connected to the ambient environment, and further analysis could illuminate key snow processes.

The PIP has been used in several field campaigns and long-term deployments all over the world, and these data are available for continued examination of snow. The PIP observations and higher-order products serve as valuable constraints on microphysical assumptions in ground-based and spaceborne remotely sensed snowfall retrievals. Additionally, this work illustrates that the higher-order PIP-derived LWE snow rate and ${\rho }_e$ values are physically consistent and valid when compared to several retrieval methods and measurements. Furthermore, the PIP captures accumulation during snow conditions under which traditional precipitation gauges are unable, due to the unique open volume design. The PIP data and higher-order products are an invaluable resource to the snow and precipitation communities.

5. Data Availability

All data used in this study can be acquired from the NASA Global Hydrology Research Center data archive (https://ghrc.nsstc.nasa.gov/home/).