What millimeter-wavelength radar reflectivity reveals about snowfall: an information-centric analysis

The ability of single-frequency, millimeter-wavelength radar reflectivity observations to provide useful constraints for retrieval of snow particle size distribution (PSD) parameters, snowfall rates, and snowfall accumulations is examined. An optimal estimation snowfall retrieval that allows analyses of retrieval uncertainties and information content is applied to observations of near-surface W-band reflectivities from multiple snowfall events during the 2006–2007 winter season in southern Ontario. Retrieved instantaneous snowfall rates generally have uncertainties greater than 100 %, but single-event and seasonal snow accumulations from the retrieval results match well with collocated measurements of accumulations. Absolute fractional differences are mainly below 30 % for individual events that have more substantial accumulations and, for the season, 12.6 %. Uncertainties in retrieved snowfall rates are driven mainly by uncertainties in the retrieved PSD parameters, followed by uncertainties in particle model parameters and, to a lesser extent, the uncertainties in the fall-speed model. Uncertainties attributable to assuming an exponential distribution are negligible. The results indicate that improvements to PSD and particle model a priori constraints provide the most impactful path forward for reducing uncertainties in retrieved snowfall rates. Information content analyses reveal that PSD slope is well-constrained by the retrieval. Given the sensitivity of PSD slope to microphysical transformations, the results show that such retrievals, when applied to radar reflectivity profiles, could provide information about microphysical transformations in the snowing column. The PSD intercept is less well-constrained by the retrieval. While applied to near-surface radar observations in this study, the retrieval is applicable as well to radar observations aloft, such as those provided by profiling ground-based, airborne, and satellite-borne radars under lighter snowfall conditions when attenuation and multiple scattering can be neglected.

1 Introduction

Radar observations focused on snowfall from platforms outside the established weather surveillance radar networks have become ubiquitous over the last 2 decades, largely due to increased interest in the role of snowfall in mid- and high-latitude microphysics, hydrology, and climate. This research accelerated with the advent of satellite-borne radars flown by missions to quantify global hydrometeor and precipitation properties. These satellite-borne radars (specifically the CloudSat mission's Cloud Profiling Radar (CPR) (Tanelli et al., 2008) and the Global Precipitation Measurement (GPM) mission's Dual-frequency Precipitation Radar (DPR) (Toyoshima et al., 2015), with two others anticipated to launch in the coming decade) are capable solely of measuring vertical profiles of radar reflectivity factor (hereafter, reflectivity) along with path-integrated attenuation under certain conditions. To understand the capabilities of these satellite-borne radars for quantifying snowfall, we must know how well radar reflectivity observations constrain snowfall properties.

To these ends, CloudSat and GPM have contributed to multiple field experiments involving ground-based radars and designed to provide, in part, ground validation data for the radar remote sensing of snowfall: the Canadian CloudSat-CALIPSO Validation Project (C3VP) (Hudak et al., 2006), the Global Precipitation Measurement (GPM) Cold-season Precipitation Experiment (GCPEx) (Skofronick-Jackson et al., 2015), the Light Precipitation Validation Experiment (LPVEx) (Petersen et al., 2011), the International Collaborative Experiment during the PyeongChang 2018 Olympics and Paralympics (ICE-POP) (Chandrasekar et al., 2019), and the Olympic Mountains Experiment (OLYMPEx) (Houze et al., 2017). These and a number of smaller, more focused field campaigns (Pettersen et al., 2020; Schirle et al., 2019) have made extensive use of small K-band profiling radars, e.g., METEK's Micro Rain Radar (MRR, Klugmann et al., 1996), but several experiments, including C3VP, GCPEx, and ICE-POP, have deployed ground-based, W-band scanning, or profiling radars. Although these ground-based radars may provide advanced capabilities such as Doppler velocity measurement, their reflectivity measurements in snowfall are a valuable resource for examining the capabilities of the satellite-borne radars (Maahn et al., 2014; Matrosov et al., 2008).

The ability of radar reflectivity to constrain snowfall properties, however, has not been well-evaluated. Snowfall exhibits a wide range of microphysical characteristics that influence radar reflectivity and snowfall rate. Most notable to casual observers are variations in particle habits: pristine dendrites, needles, columns, plates, and bullets; aggregates of the same; pellets and graupel for example. Underlying these differences in habit are variations in mass, and, given a particular mass, variations in how mass is distributed within the particle. Unlike longer-wavelength radars for which radar backscattering properties of snow particles are sensitive primarily to particle mass, at millimeter wavelengths those properties are additionally sensitive to particle shape. Investigations of particle mass and area (an aspect of shape) (Kajikawa, 1972, 1975, 1982; Zikmunda and Vali, 1972, 1977; Heymsfield, 1972; Locatelli and Hobbs, 1974; Mitchell et al., 1990; Mitchell, 1996; Heymsfield and Miloshevich, 2003) have painstakingly determined the broad extent of these variations. Along with differences in single-particle properties, populations of falling snow particles vary substantially in their concentrations with size (i.e., the spectral particle size distribution, PSD) based on measurements from the ground (e.g., Nakada and Terada, 1935; Imai et al., 1955; Gunn and Marshall, 1958; Rogers, 1973; Brandes et al., 2007) and, more recently, with the advent of imaging particle probes, from aircraft (e.g., Passarelli, 1978; Gordon and Marwitz, 1984, 1986; Braham, 1990; Woods et al., 2008; Heymsfield et al., 2008, and references therein). The observed particle concentrations vary over several orders of magnitude.

In radar-based remote sensing scenarios when these properties are not known, these variations produce uncertainty in the relationship between radar reflectivity factor (hereafter, reflectivity) and associated water content and snowfall rate. A common approach to estimating this uncertainty has been to evaluate modeled reflectivities, water contents, and snowfall rates using a range of assumed particle models and PSDs. The results are often expressed using relationships between reflectivity and snowfall rate (“Z–S” relationships) (Liu, 2008; Kulie and Bennartz, 2009; Matrosov et al., 2008, 2009). This approach allows the uncertainty in a retrieved snowfall rate to be estimated, but the existing studies have not provided insight into the dominant sources of uncertainty nor into the ability of observed radar reflectivity to constrain various properties controlling the snowfall rate. Posselt et al. (2015) examined uncertainties and information content for radar observations of mixed- and ice-phase regions of a convective storm but targeted radar systems with more advanced capabilities. Mascio and Mace (2017) used CloudSat and aircraft observations to assess how uncertainties in the ice particle mass-dimension relationship contribute to radar reflectivity forward model uncertainties but used known, observed particle size distributions and did not examine the influence of the mass-dimension uncertainties on snowfall retrieval performance.

In this work, we provide uncertainty and information content analyses for retrieving snowfall from observations of radar reflectivity at millimeter wavelengths, focusing on W-band (94 GHz). The results are representative of the general problem of estimating snowfall from such remote radar reflectivity observations without supplementary collocated observations of snow particle mass-dimension relationships, fall speeds, and particle size distributions. The results apply particularly to observations by the CPR (Tanelli et al., 2008) and by the DPR's Ka-band radar, but also to reflectivity measurements from ground-based radars such as the MRR (Klugmann et al., 1996) and the Department of Energy Atmospheric Radiation Measurement (ARM) program's Millimeter Wavelength Cloud Radar (Moran et al., 1998) and Ka-band ARM Zenith Radar (KAZR) (Bharadwaj et al., 2013). The retrieval method used here is the foundation for the retrieval used for CloudSat's 2C-SNOW-PROFILE product; that application is the subject of a future companion paper. Our objectives here are to identify the snowfall properties that are best constrained by such observations and the most significant sources of uncertainty in the radar retrieval of snowfall. The results establish a performance baseline for reflectivity-only observations of snowfall, indicate where uncertainty reduction efforts should be focused, and suggest what improvements to radar-observing systems could be most beneficial.

The analyses use the optimal estimation (OE) retrieval technique (Rodgers, 2000), which inherently diagnoses information content and uncertainties in retrieved quantities subject to specified uncertainties in measurements, forward models, and a priori knowledge of the quantities to be retrieved (L'Ecuyer et al., 2006; Cooper et al., 2006). The retrieval produces best estimates of snow size distribution parameters by using the radar reflectivity observations to refine a priori estimates of those parameters (Sect. 2). The information content metrics provided by OE require all sources of uncertainties in the retrieval process to be specified. These are discussed in Sect. 3. Ground-based radar and precipitation observations allow the retrieval to be tested, showing that size distribution width is best constrained by the retrieval and that uncertainties in retrieved size distribution parameters (but not uncertainties due to the assumed exponential form of the PSD itself) are the strongest contributors to uncertainties in estimated snowfall rates (Sect. 4). The results suggest that the retrieved size distribution widths could be useful for diagnosing changes in PSD resulting from microphysical processes (Lo and Passarelli, 1982) and that improved observational constraints on size distribution parameters, as might be provided by dual-wavelength radar observations (Matrosov, 1998), would likely enhance snowfall retrieval performance (Sect. 5).

2 Retrieval method

The retrieval uses measurements of reflectivity to estimate snow microphysical properties and to quantify water content and snowfall rate. At the wavelengths characteristic of cloud radars such as CloudSat and shorter-wavelength precipitation radars, scattering by precipitation-sized particles does not follow the Rayleigh approximation, and both attenuation and multiple scattering may affect the radar signal. At these wavelengths, snow particle scattering and extinction properties depend not only on mass, but on shape as well. With even simple parameterized expressions for particle mass, shape, and size distribution, single-frequency observations of radar reflectivity alone are insufficient to reasonably constrain the resulting set of parameters.

To address this insufficiency, retrievals must incorporate a priori information about particle microphysical and scattering properties. This is accomplished here using OE (Rodgers, 2000), a Bayesian technique that allows a priori information to be included explicitly. The input for this retrieval is the Z_e observed by the radar for a range gate identified as containing snow. For notational consistency with other work, we show this as a vector:

$$\begin{array}{}\text{(1)}& \mathit{y}=\left[Z_{e,\mathrm{1}}^{\mathrm{obs}}\right].\end{array}$$

A forward model $\mathit{F}(\mathit{x},\stackrel{\mathrm{̃}}{\mathit{b}})$ relates y to x, a state vector of unknown properties to be retrieved, as

$$\begin{array}{}\text{(2)}& \mathit{y}=\mathit{F}(\mathit{x},\stackrel{\mathrm{̃}}{\mathit{b}})+\mathit{ϵ},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathit{b}}$ are parameters not being retrieved but which influence the forward model results. The forward model approximates the true physical relation between x and y, and there are uncertainties associated with both the observations y and the forward model parameters $\stackrel{\mathrm{̃}}{\mathit{b}}$. ϵ represents the total uncertainty due to all sources. OE attempts to find $\widehat{\mathit{x}}$, an estimate of the state which maximizes the posterior conditional probability density function (PDF) P(x|y), subject also to prior knowledge about the values of x. This prior knowledge is described by expected values x_a and their covariances S_a. Assuming Gaussian statistics for the model-measurement errors and the a priori state, minimizing the cost function

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Phi }(\mathit{x},\mathit{y},{\mathit{x}}_a)& ={\left(\mathit{y}-\mathit{F}(\mathit{x},\stackrel{\mathrm{̃}}{\mathit{b}})\right)}^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\left(\mathit{y}-\mathit{F}(\mathit{x},\stackrel{\mathrm{̃}}{\mathit{b}})\right)\\ & +{\left(\mathit{x}-{\mathit{x}}_a\right)}^{\mathsf{T}}{\mathbf{S}}_a^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_a\right)\end{array}\end{array}$$

with respect to x gives this PDF, where S_ϵ is the covariance matrix representing the uncertainties ϵ. The Gaussian assumption is reasonable if the expected values and covariance matrices are known for the model-measurement uncertainties and the a priori state, but other details are lacking. In that case, the Gaussian form maximizes the entropy of a PDF (Shannon and Weaver, 1949; Rodgers, 2000). Assuming an alternate form would introduce constraints on the retrieval that are not justified based on the limited knowledge of the PDF.

Provided the forward model is not excessively nonlinear, Newtonian iteration

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\widehat{\mathit{x}}}_{i+\mathrm{1}}& ={\widehat{\mathit{x}}}_i+{\left({\mathbf{S}}_a^{-\mathrm{1}}+{\mathbf{K}}_i^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i\right)}^{-\mathrm{1}}\\ & \left[{\mathbf{K}}_i^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\left(\mathit{y}-\mathit{F}({\widehat{\mathit{x}}}_i,\stackrel{\mathrm{̃}}{\mathit{b}})\right)-{\mathbf{S}}_a^{-\mathrm{1}}({\widehat{\mathit{x}}}_i-{\mathit{x}}_a)\right]\end{array}\end{array}$$

leads to $\widehat{\mathit{x}}$, where K is the Jacobian of the forward model with respect to x, and ${\mathbf{K}}_i=\mathbf{K}\left({\widehat{\mathit{x}}}_i\right)$. Iteration continues until the squared difference in successive ${\widehat{\mathit{x}}}_i$ normalized by the current estimate of the a posteriori covariance ${\widehat{\mathbf{S}}}_x$ is much smaller than the number of state vector elements (Rodgers, 2000). At convergence, this covariance of $\widehat{\mathit{x}}$ is obtained as

$$\begin{array}{}\text{(5)}& {\widehat{\mathbf{S}}}_x={\left({\widehat{\mathbf{K}}}^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\widehat{\mathbf{K}}+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}},\end{array}$$

where $\widehat{\mathbf{K}}=\mathbf{K}\left(\widehat{\mathit{x}}\right)$. As a diagnostic test of the results, a χ² statistic is calculated using the retrieved state vector in Eq. (3). A value near the number of observations suggests correct convergence (Marks and Rodgers, 1993). Several metrics, determined from the retrieved state and based on information theory, provide insight into the retrieval performance; these metrics are presented in Sect. 4.

2.1 The forward model

To assess the information provided purely by reflectivity observations, whether from ground-, aircraft-, or space-based radars, the retrieval ignores attenuation and multiple scattering. This treatment would be appropriate for cases with little intervening scattering and extinction between the radar and observed snowfall, such as when the radar bin containing the snowfall of interest is near the radar or under light snowfall conditions. For such a case, the singly-scattered reflectivity $Z_e^{\mathrm{ss}}$ as a function of range R from the radar is given by

$$\begin{array}{}\text{(6)}& Z_e^{\mathrm{ss}}\left(R\right)={\displaystyle \frac{{\mathrm{\Lambda }}^{\mathrm{4}}}{{∥K_w∥}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{5}}}}\underset{D_{\min }}{\overset{D_{\max }}{\int }}N(D,R){\mathit{\sigma }}_{\mathrm{bk}}(D,R)\mathrm{d}D,\end{array}$$

where σ_bk(D,R) is the backscatter cross section for particle size D at range R, N(D,R) is the PSD at range R, Λ is the radar wavelength, and K_w is the dielectric factor for water.

2.1.1 Forward model parameters: snow particle model

Backscattering and extinction cross sections depend intimately on particle mass, shape, and orientation relative to the radar beam. These properties are highly variable for snow particles, and the approach used here is to specify their PDFs a priori using best estimates and treat their variability as a source of uncertainty in the retrieval. We adopt the common model (e.g., Locatelli and Hobbs, 1974; Mitchell, 1996) in which mass and horizontally projected area are described using power laws,

$$\begin{array}{}\text{(7)}& {\displaystyle }m\left(D_{\mathrm{M}}\right)=\mathit{\alpha }{D}_{\mathrm{M}}^{\mathit{\beta }},\text{(8)}& {\displaystyle }{A}_{\mathrm{p}}\left(D_{\mathrm{M}}\right)=\mathit{\gamma }{D}_{\mathrm{M}}^{\mathit{\sigma }},\end{array}$$

on particle maximum dimension, D_M, and use the particle properties and shape “B8pr-30” (Wood et al., 2015), an idealized branched spatial particle that was found to minimize bias in simulated reflectivities versus coincident W-band radar observations. That work used in situ measurements and remotely sensed X-band reflectivity observations of snow from C3VP (Hudak et al., 2006) along with previously reported single-particle measurements to develop best estimates and covariances for the power-law parameters α, β, γ, and σ. These results then constrained discrete dipole approximation calculations using DDSCAT (Draine and Flatau, 1994) to obtain best estimates of snow particle single-scattering properties and their uncertainties at the desired wavelengths. These a priori descriptions of size-resolved particle mass, A_p, σ_bk, σ_ext, and their uncertainties constitute the particle model used in the retrieval and are summarized in Appendix B.

2.2 The retrieved state

The relationship described by Eq. (6) requires information about particle size distributions and single-scattering properties. With scattering properties and their uncertainties specified a priori as described in Sect. 2.1.1, this leaves the snow PSD parameters and their PDFs to be determined by the retrieval.

Snow PSDs are frequently characterized as exponential:

$$\begin{array}{}\text{(9)}& N\left(D\right)=N_{\mathrm{0}}\exp \left(-\mathit{\lambda }D\right),\end{array}$$

where λ is the slope of the distribution and N₀ its intercept. Rogers (1973) used photographs of snowflakes to develop estimates of snow size distributions based on actual dimensions and found snow size distributions to be exponential. Brandes et al. (2007) evaluated both exponential and gamma forms, which have the ability to represent sub- or super-exponential behavior, for snow size distributions observed by a 2D video disdrometer over the course of several winter seasons. Although about 22 % of the observed snow distributions exhibited super-exponential features, more commonly the fitted gamma distributions were nearly equivalent to exponential distributions. Several aircraft-based studies using in situ observations under a wide range of atmospheric conditions have confirmed exponential behavior, especially at larger particle sizes (Passarelli, 1978; Houze et al., 1979; Lo and Passarelli, 1982; Gordon and Marwitz, 1984; Braham, 1990; Woods et al., 2008). While other studies of aircraft observations have noted departures from exponential behavior (e.g., “super-” or “sub-exponential”, Herzegh and Hobbs, 1985), Heymsfield et al. (2008) examined the suitability of exponential distributions for snow. They found that fitted exponential distributions, when used to simulate IWC and Ze, could provide generally good agreement with IWC and Ze calculated directly from the observed discrete size distributions. These studies support the adequacy of exponential distributions for retrieving snowfall. D may be an actual dimension of the snow particle, the diameter of an equivalent mass ice sphere, or the melted drop diameter. The choice is significant because N₀ and λ depend on the choice of D. For this work, we use the maximum particle dimension, D_M, because D_M is closely related to the particle dimensions measured by imagers such as video disdrometers (Wood et al., 2013) and aircraft particle probes, making comparisons with other datasets more straightforward.

The exponential size distribution parameters N₀ and λ are the desired state variables. Values for N₀ may range over several orders of magnitude, so log (N₀) is retrieved instead. The variability of λ is significantly smaller than that of N₀; however, examination of fitted exponential distributions from C3VP snow events indicated that the distribution of values for λ was strongly non-Gaussian. The log-transformed values are much less skewed (Fig. 1a), and accordingly, log (λ) is retrieved instead. The corresponding state vector to be retrieved is then

$$\begin{array}{}\text{(10)}& \widehat{\mathit{x}}=\left[\begin{array}{l}\log \left(N_{\mathrm{0}}\right)\\ \log \left(\mathit{\lambda }\right)\end{array}\right],\end{array}$$

and the associated covariance matrix obtained from the retrieval is of the form

$$\begin{array}{}\text{(11)}& {\widehat{\mathbf{S}}}_x=\left[\begin{array}{ll}{s}^{\mathrm{2}}\left(\log \left(N_{\mathrm{0}}\right)\right)& \begin{array}{c}s\left(\log \left(N_{\mathrm{0}}\right),\right.\\ \left.\log \left(\mathit{\lambda }\right)\right)\end{array}\\ \begin{array}{c}s\left(\log \left(N_{\mathrm{0}}\right),\right.\\ \left.\log \left(\mathit{\lambda }\right)\right)\end{array}& s^{\mathrm{2}}\left(\log \left(\mathit{\lambda }\right)\right)\end{array}\right].\end{array}$$

Figure 1 (a) Histograms of λ and log (λ) fitted to C3VP SVI observations. (b) Estimates of λ and N₀ determined from fits to size distributions from C3VP observations, with values provided from several earlier studies for comparison.

2.3 Prior estimates of the state

For each profile, the a priori state consists of a vector of expected values x_a and the corresponding covariance matrix S_a, having the same sizes as the state vector x (Eq. 10) and its covariance matrix S_x (Eq. 11). A priori estimates of log (N₀) and log (λ) are determined using temperature-based parameterizations derived using snow PSDs observed during C3VP and other field experiments. Exponential size distributions were fit to the observed size spectra from both ground-based Snowflake Video Imager, or SVI (Newman et al., 2009; Wood et al., 2013), and from 2D particle probes carried aboard the National Research Council Canada's Convair-580 during three C3VP research flights (Fig. 1b). Results from a number of earlier studies are shown as well for comparison, including ground-based observations taken in and near the Rocky Mountain Front Range (Rogers, 1973; Brandes et al., 2007); and aircraft observations over the central Sierra Nevada (Gordon and Marwitz, 1984, 1986), in lake effect snow over Lake Michigan (Braham, 1990), in synoptic snowfall over central Illinois (Passarelli, 1978), and in orographic and frontal wintertime precipitation in the Pacific Northwest (Woods et al., 2008). Also shown are similar fits performed on 2D probe observations from a Wakasa Bay research flight on 27 January 2003 (Lobl et al., 2007). The results suggest that the C3VP observations adequately represent snowfall from a number of different regimes, although the number concentrations from several studies are at the margins of the C3VP observations.

Both λ and N₀ have been observed to vary log-linearly with temperature (e.g., Houze et al., 1979; Woods et al., 2008; and works reviewed in Ryan, 1996). Fits were therefore constructed for both parameters using the combined C3VP aircraft and SVI data and uncertainties estimated using residual standard deviations (RSDs) calculated for data binned into 2 K intervals (Fig. 2). The narrow temperature ranges for the Wakasa Bay and Brandes et al. (2007) observations make comparisons against the C3VP temperature dependence uninformative. For λ, the Rogers (1973) observations are largely outside the bounds of the RSDs but are generally consistent with the C3VP histogram for warmer temperatures. The aircraft observations other than Wakasa Bay follow a temperature trend similar to the C3VP observations. For N₀, several of the comparison datasets lie mostly above the RSD bounds but would be well within a +2 RSD bound.

Figure 2 Dependence of log (λ) and log (N₀) on temperature. Central red lines show the best-fit relationships, while the upper and lower blue lines show bounds given by ±1 residual standard deviation. The shaded gray shows the 2D histogram of values for the C3VP surface and aircraft observations (a and c). Symbols (b and d) match those from Fig. 1 except that, in lieu of symbols for Woods et al. (2008), the dashed black line shows a linear best fit reported by the authors.

Based on the similarity of C3VP to results from other experiments, the a priori states derived from these observations can be expected to represent a broad range of snowfall regimes and were adopted for the retrieval. A priori values for log (λ) and log (N₀) were estimated from the linear fits as

$$\begin{array}{}\text{(12)}& \begin{array}{l}\log \left({\mathit{\lambda }}_{\mathrm{ap}}\right)=-\mathrm{0.03053}(T-\mathrm{273}.)-\mathrm{0.08258},\\ \log \left(N_{\mathrm{0},\mathrm{ap}}\right)=-\mathrm{0.07193}(T-\mathrm{273}.)+\mathrm{2.665},\end{array}\end{array}$$

with λ in mm⁻¹, N₀ in m⁻³ mm⁻1, and T in K. The RSDs show little variation with temperature except in the vicinity of 240 K, where they increase substantially. These large RSDs are in response to a few outlying samples with small λ and N₀ values. Accordingly, variances were treated as constant and were estimated as the squared RSDs averaged over all temperatures. The uncertainty model also includes the covariance between log (N₀) and log (λ). Correlation coefficients were evaluated for each of the temperature-binned data subsets, giving a mean coefficient of 0.72 with a standard deviation of 0.12. The a priori covariance was modeled as $\mathrm{0.72}\cdot s\left(\log \left({\mathit{\lambda }}_{\mathrm{ap}}\right)\right)\cdot s\left(\log \left(N_{\mathrm{0},ap}\right)\right):$

$$\begin{array}{}\text{(13)}& \begin{array}{l}{s}^{\mathrm{2}}\left(\log \left({\mathit{\lambda }}_{\mathrm{ap}}\right)\right)=\mathrm{0.133},\\ s^{\mathrm{2}}\left(\log \left(N_{\mathrm{0},\mathrm{ap}}\right)\right)=\mathrm{0.95},\\ s\left(\log \left({\mathit{\lambda }}_{\mathrm{ap}}\right),\log \left(N_{\mathrm{0},ap}\right)\right)=\mathrm{0.26}.\end{array}\end{array}$$

3 Implementation and uncertainty sources

Applying the exponential distribution in Eq. (6), the singly-scattered non-attenuated reflectivity $Z{}_e^{\mathrm{ss}}$ is

$$\begin{array}{}\text{(14)}& \begin{array}{rl}Z{}_e^{\mathrm{ss}}\left(R\right)& ={\displaystyle \frac{{\mathrm{\Lambda }}^{\mathrm{4}}}{{∥K_w∥}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{5}}}}\underset{D_{\mathrm{M},\min }}{\overset{D_{\mathrm{M},\max }}{\int }}{N}_{\mathrm{0}}\exp \left(-\mathit{\lambda }{D}_{\mathrm{M}}\right){\mathit{\sigma }}_{\mathrm{bk}}\\ & (D_{\mathrm{M}},\stackrel{\mathrm{̃}}{\mathit{b}})\mathrm{d}{D}_{\mathrm{M}}.\end{array}\end{array}$$

The backscatter cross section σ_bk has been written to show its dependence on a vector of parameters $\stackrel{\mathrm{̃}}{\mathit{b}}$ as well as on D_M. The vector $\stackrel{\mathrm{̃}}{\mathit{b}}$ includes the parameters for the mass- and area-dimension relations α, β, γ, and σ which were used to construct the particle models from which the scattering properties were calculated. The tilde indicates that these parameters are approximations of the true values and a source of uncertainty.

3.1 Model-measurement uncertainties

The error covariance matrix S_ϵ is

$$\begin{array}{}\text{(15)}& \begin{array}{ll}{\mathbf{S}}_{\mathit{ϵ}}& ={\mathbf{S}}_y+{\mathbf{S}}_F\\ & ={\mathbf{S}}_y+{\mathbf{S}}_B^{\mathrm{ss}}+{\mathbf{S}}_F^{\mathrm{ss}},\end{array}\end{array}$$

where S_y is the covariance matrix for the measurement uncertainties and S_F is that for the singly-scattered reflectivities given in Eq. (14). The forward-model uncertainties may be further decomposed as the sum of two terms: ${\mathbf{S}}_B^{\mathrm{ss}}$, which is a covariance matrix describing uncertainties due to the forward model parameters $\stackrel{\mathrm{̃}}{\mathit{b}}$, and ${\mathbf{S}}_F^{\mathrm{ss}}$, which is a covariance matrix describing uncertainties due to other assumptions in the calculation of $Z_e^{\mathrm{ss}}$.

3.1.1 Uncertainties for measured reflectivities

The sources of reflectivity measurement error include errors in the absolute radiometric calibration and measurement noise. For this work, we assume the radar is well-calibrated, leaving noise as the uncertainty source. To estimate S_y, we model the noise using the well-characterized CloudSat CPR (Tanelli et al., 2008). For reflectivities above −10 dBZ, 1 standard deviation of noise as a fraction of the mean signal is about −16 dB, while for reflectivities below −10 dBZ, noise is an increasing fraction of the signal, reaching 0 dB at the minimum detectable signal of −30 dBZ (Richard T. Austin, personal communication, 4 November 2008). The resulting uncertainties range from 3 dBZ for a reflectivity of −30 dBZ to about 0.1 dBZ for reflectivities above −10 dBZ_e (Fig. 3).

Figure 3 Estimated measurement uncertainty, based on 1 standard deviation of noise for the CloudSat CPR.

3.1.2 Forward model uncertainties

Uncertainties ${\mathbf{S}}_B^{\mathrm{ss}}$ due to the forward model parameters $\stackrel{\mathrm{̃}}{\mathit{b}}={\left(\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\sigma }\right)}^{\mathtt{T}}$ that describe the snow particle model were examined in Wood et al. (2015) as

$$\begin{array}{}\text{(16)}& {\mathbf{S}}_B^{\mathrm{ss}}={\mathbf{K}}_b{\mathbf{S}}_b{\mathbf{K}}_b^{\mathsf{T}},\end{array}$$

where K_b is the Jacobian of the forward model reflectivities with respect to the parameters $\stackrel{\mathrm{̃}}{\mathit{b}}$ and S_b is the covariance matrix for the parameters. The Jacobian K_b depends on the estimated state $\widehat{\mathit{x}}$ and is computed at each iterative step. At each step, the forward model is used to calculate reflectivity perturbations that result from perturbations of the parameters α, β, γ, and σ. The ratio of each reflectivity perturbation to its parameter perturbation gives an element of the Jacobian. The parameter perturbation affects the reflectivity via changes to the corresponding particle scattering properties. The perturbed scattering properties are precomputed with DDSCAT (Draine and Flatau, 1994) by using the perturbed parameter to generate discrete dipole models following the process described in Wood et al. (2015). Wood et al. (2015) found the resulting forward model uncertainties to be near 5 dB, increasing to as high as 15 dB for very broad distributions.

${\mathbf{S}}_F^{\mathrm{ss}}$ quantifies uncertainties due to other assumptions and limitations in the forward model reflectivity calculation. Wood et al. (2015) looked at uncertainties due to the random component of dipole placement within discrete dipole approximation (DDA) models for a particular particle shape and found them negligible. Other sources include the assumption of the shape of the distribution as exponential, the choice of particle shape, and the discretization and truncation of the integrations over size distribution.

Errors due to the assumed exponential shape were evaluated using a dataset of 4080 SVI-measured, discrete, 5 min-long snow PSDs from C3VP. Simulated reflectivities and snowfall rates were calculated using the B8pr-30 particle model and the Mitchell and Heymsfield (2005) terminal velocity model. Exponential distributions were fit to the observed discrete PSDs using orthogonal distance regression (Boggs et al., 1992; Jones et al., 2001) with uncertainty estimates per Wood et al. (2013). The fitted distributions were scaled in number concentration to match the snowfall rates simulated from the discrete distributions. The fitted distributions were then used to simulate reflectivities for comparison against those from the discrete distributions. Errors are negligible at high reflectivities but increase as reflectivity decreases (Fig. 4). Bias is negligible, and the total uncertainty is modeled as

$$\begin{array}{}\text{(17)}& s_F^{\mathrm{2}}\left(\mathrm{dB}\right)={\left[\exp \left(-\left({\mathrm{dBZ}}_e+\mathrm{14}\right)/\mathrm{16}\right)\right]}^{\mathrm{2}},\end{array}$$

reaching a maximum of 1 dB of uncertainty at −15 dBZ_e.

Figure 4 Actual rms errors and the fitted model for uncertainty due to the assumed exponential size distribution.

Uncertainties due to shape were evaluated using the same SVI dataset to which the alternate particle shapes Ep (ellipsoidal) and B8pr-45 (branched spatial particle with a larger aspect ratio than B8pr-30) from Wood et al. (2015) were applied to simulate reflectivities. These alternate shapes are constrained to have the same mass-dimension relationship as used for the B8pr-30 particle model used in this work, so differences are due only to particle shape. Figure 5 shows total and variance-only rms errors. From these results we estimate the shape uncertainty to be 2 dB.

Figure 5 Errors in reflectivity for the Ep and B8pr-45 shapes compared to the B8pr-30 shape. Errors shown are total (bias + variance) and variance only.

Finally, truncation and discretization errors were evaluated using the same SVI PSD dataset. These are errors that result from the discrete treatment of the integrations over size distribution, errors due to both the limited maximum D_M in the particle model and in the limited resolution of the particle model. Truncation errors were evaluated using analytic exponential PSDs fitted to the SVI PSD dataset as described previously. The particle model backscatter properties were augmented to D_M=40 mm by linearly extrapolating backscatter efficiencies, and then reflectivities were calculated using integrations to both the standard (maximum D_M=18 mm) and augmented size ranges. The bias and scatter of the truncation errors were −0.1 and 0.42 dB. To evaluate discretization errors, a high-resolution version of the particle model backscatter properties was created by interpolating backscatter efficiencies so that the particle size resolution of the particle model was increased by a factor of 2. Reflectivities were then calculated and compared against those from the standard-resolution particle model. The bias and scatter of the discretization errors were 0.00 and 0.02 dB.

Figure 6 Histograms of errors for truncation and discretization.

3.2 Snowfall rate and uncertainties

The snowfall rate P in units of liquid water depth per unit time is

$$\begin{array}{}\text{(18)}& P\left(R\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{liq}}}}\underset{D_{\mathrm{M},\min }}{\overset{D_{\mathrm{M},\max }}{\int }}N(D_{\mathrm{M}},R)m\left(D_{\mathrm{M}},R\right)V\left(D_{\mathrm{M}},R\right)\mathrm{d}{D}_{\mathrm{M}},\end{array}$$

where m(D_M,R) is particle mass, V(D_M,R) is fall speed, and ρ_liq is the density of liquid water. Particle mass is provided by Eq. (7). Fall speed is assumed to equal terminal velocity, which is calculated from the model of Mitchell and Heymsfield (2005) using particle mass, the horizontally projected area from Eq. (8), and environmental pressure and temperature from collocated observations. Uncertainties for the estimated snowfall rate are determined in a manner similar to that used for the forward model uncertainties. The total variance S_P is decomposed as

$$\begin{array}{}\text{(19)}& {\mathbf{S}}_P={\mathbf{S}}_{\widehat{x},P}+{\mathbf{S}}_{\stackrel{\mathrm{̃}}{b},P}+{\mathbf{S}}_{v,P}+{\mathbf{S}}_{\exp ,P},\end{array}$$

where the terms on the right represent the variances resulting from (1) retrieved state uncertainties, (2) particle model parameter uncertainties, (3) uncertainties in the fall-speed model and its parameters, and (4) assuming an exponential form for the PSD, respectively.

Contributions from uncertainties in the retrieved state and in the particle model parameters are determined using linearized error propagation (e.g., following a form like Eq. 16). For ${\mathbf{S}}_{\stackrel{\mathrm{̃}}{b},P}$, which gives the snowfall rate variance that results from uncertainties in the particle model parameters, this means that the Jacobian ${\mathbf{K}}_{\stackrel{\mathrm{̃}}{b},P}$ is calculated for the snowfall rate with respect to the particle model parameters α, β, γ, and σ, following the process described for the reflectivity Jacobian in Sect. 3.1.2. Then

$$\begin{array}{}\text{(20)}& {\mathbf{S}}_{\stackrel{\mathrm{̃}}{b},P}={\mathbf{K}}_{\stackrel{\mathrm{̃}}{b},P}{\mathbf{S}}_b{\mathbf{K}}_{\stackrel{\mathrm{̃}}{b},P}^{\mathsf{T}},\end{array}$$

where S_b is the covariance matrix for the particle model parameters as determined in Wood et al. (2015).

Fall-speed contributions are handled following Wood et al. (2014). Snowfall rate uncertainties due to the assumed exponential form of the size distribution are determined using the SVI PSD dataset in an approach analogous to that for Eq. (17). In this approach, number concentrations for the fitted exponential distributions were scaled so that reflectivities were matched, and then snowfall rate errors were evaluated. The fractional uncertainty in snowfall rate was found to be

$$\begin{array}{}\text{(21)}& f_{\mathrm{P}}=-\mathrm{0.06}\log \left(P\right)+\mathrm{0.05},\end{array}$$

from which the necessary variance can be determined. Uncertainties from each of the four sources are treated as uncorrelated.

4 Retrieval performance tests with ground-based radar observations

During C3VP, a vertically pointing W-band radar (the Jet Propulsion Laboratory's Airborne Cloud Radar, ACR) was deployed on the ground at CARE. In all, about 28 h of ACR radar profiles of snowfall were recorded at approximately 2.8 s intervals. These observations represent 17 distinct snow events that occurred over 18 d between 3 November 2006 and 2 March 2007; however, most of the accumulations were concentrated during nine of the events (Table 1). These observations include portions of three of the cases that were used to develop the snow particle microphysical models (cases SYN1, LES1, and LES2, Wood et al., 2015). Of the nearly 36 000 ACR profiles in these observations, approximately 7300 are from cases SYN1, LES1, and LES2. Further, as described in Wood et al. (2015), ACR reflectivities from 12 of the events between 2 December 2006 and 26 February 2007 were used to constrain the snow particle models' scattering properties to give unbiased reflectivities. This overlap should be kept in mind when evaluating the retrieved snowfall rates and estimated accumulation, but it should not substantially affect the assessments of retrieval uncertainties, uncertainty sources, and information content metrics that follow.

The retrieval was applied to the ACR reflectivities observed in the single range bin nearest the surface, at 197 m above ground level (AGL). Temperatures and pressures needed by the retrieval to perform snow detection, calculate fall speeds, and establish the a priori states were obtained from nearby surface meteorology observations. Because of the short distance to the target range bin, attenuation along the path was neglected. The retrieved snowfall rates produce a Z–S relationship that is most similar to that developed by Kulie and Bennartz (2009) for an aggregate particle model denoted as the Hong (2007) aggregate (HA) (Fig. 7). For warmer temperatures and mid-range reflectivities, the Z–S relationship becomes more similar to that of Liu (2008) and the LR3 relationship of Kulie and Bennartz.

Figure 7 Z–S values as a function of temperature for this retrieval compared against those from M07, Matrosov (2007); L08, Liu (2008); and KB09_LR3, KB09_HA, and KB09_SS, Kulie and Bennartz (2009).

For comparisons, snowfall rate observations were obtained at 1 min intervals from the Vaisala FD12P (Vaisala Oyj, 2002) and scaled to provide unbiased accumulations relative to the nearby Dual Fence Intercomparison Reference, or DFIR (Goodison et al., 1998). The retrieved ACR snowfall rates, P_ACR, were matched to the nearest-in-time observed snowfall rate, P_FD12P.

Table 1 Accumulations by event for the ACR retrievals. Duration shows the elapsed time of ACR observations for which retrievals were performed. Fractional differences are relative to FD12P accumulations.

Two distinct events, indicated as ^a and ^b, occurred on 20 January 2007. ^c Accumulations adjusted to remove anomalies indicated in Fig. 8.

Time series of P_ACR and P_FD12P show a high degree of agreement over most of the observing period (Fig. 8a). This is not extraordinary given the dependence of the retrieval's particle microphysical and scattering properties on portions of the C3VP data. Two notable exceptions occur near time indices 25 000 and 32 500, however, when the FD12P recorded snowfall rates above 1 mm LWE h⁻¹, while the retrieved values are substantially smaller. Examining the time series of ACR reflectivities shows that the ACR did not observe high reflectivities during these periods (Fig. 8b). The first of these anomalies occurred on 22 February 2007 from 11:20 to 12:05 UTC, while the second occurred on 1 March 2007 between 22:15 and 22:50 UTC. For both, the ACR operator made note of the heavy snowfall, suggesting that both the FD12P and the ACR observed similar snowfall rates. Based on soundings, Environment Canada forecasts, and ACR operator observations, these anomalies appear to correspond to melting aloft, ice pellets, and freezing rain (Wood, 2011). These conditions could also have been favorable for formation of large, heavy aggregates. It seems likely that the conditions produced snowfall whose properties were strongly inconsistent with the particle properties assumed in the retrieval, although reflectivities did not change substantially.

Figure 8 (a) Time series of snowfall rates retrieved from ACR reflectivities and observed. (b) Corresponding time series of ACR reflectivities. Each time index indicates a 2.8 s observation by the ACR. Snowfall rates retrieved for the ACR used the reflectivity in the range bin nearest the surface, at 197 m AGL.

Accumulations were calculated from both P_ACR and P_FD12P with and without the two anomalies described above (Fig. 9). Accumulations agree substantially during the first 16 h but diverge somewhat beyond that, again noting the dependence of the retrieval's assumed microphysical and scattering properties on portions of the C3VP data. With the anomalies included the final difference between the accumulations is 2 mm. With the anomalies removed that difference is reduced to 0.7 mm. For individual events, absolute fractional differences between the ACR and FD12P accumulations can range to 50 % and upwards (Table 1), but these large values are associated mainly with events with small accumulations. For events with larger accumulations, the absolute fractional differences are mostly below 30 %. At seasonal timescales, the random components in event-total accumulations are likely uncorrelated, leading to offsetting errors when calculating seasonal accumulations. The time series of absolute fractional differences between the ACR-derived and FD12P accumulations begins with large fractional differences. Within 5 h and over the initial three events, the fractional differences reduce to less than 5 % and then remain below 20 % for the remainder of the season.

Figure 9 Snow accumulations computed from P_ACR and P_FD12P. The accumulations are for 17 snow events observed by the ACR on 18 d between 3 November 2006 and 2 March 2007, but accumulations are principally from nine events (Table 1). The events were concatenated sequentially in time, and the time axis indicates the cumulative time over all events. (a) Accumulations from all observations and corresponding retrieval results, (b) accumulations with two anomalous periods identified in Fig. 8 removed, and (c) fractional differences in accumulations shown in (b), with distinct colors indicating individual events.

4.1 Snowfall rate uncertainties

Uncertainties in instantaneous snowfall rate estimates, taken to be the square root of the total variance evaluated as shown in Eq. (19), were evaluated by binning the fractional uncertainties by snowfall rate and then averaging and taking standard deviations. Mean fractional uncertainties range from 150 % to 185 %, and the range for ±1 standard deviation extends from about 145 % to 190 % (Fig. 10). The fractional uncertainties generally increase with increasing snowfall rate, but above 0.5 mm LWE h⁻¹ the means and standard deviations diminish and result from only a small number of samples in each bin. For comparison, uncertainties for FD12P precipitation rates at 5 min resolution were estimated at 0.03 mm h⁻¹ for rates less than 0.05 mm h⁻¹, 50 % for rates up to 0.5 mm h⁻¹, and 30 % for rates larger than 0.5 mm h⁻¹ by Wood et al. (2014) based on comparisons against the Precipitation Occurrence Sensor System.

Figure 10 Instantaneous fractional uncertainties in snowfall rate. The central line shows mean fractional uncertainties and the error bars show ±1 standard deviation.

To evaluate the importance of each source of uncertainty, variances from each of the sources from Eq. (19) (retrieved state, microphysical parameters, fall-speed parameterization, or exponential distribution) were extracted separately, and then fractions of total variance were calculated. To allow the trends in each source to be shown as a function of snowfall rate (Fig. 11), the fractions were binned by snowfall rate and averaged. As snowfall rates increase up to 0.5 mm h⁻¹, the variance due to the retrieved state becomes a more significant contributor to the total variance, while the contributions from the other sources diminish. The contribution due to the assumed exponential PSD shape is not significant.

Figure 11 Instantaneous fractional variances for snowfall rate resolved by source.

The instantaneous uncertainties for snowfall rate include uncertainties due to random errors and biases in the retrieval components and observations. For accumulations or mean rates evaluated over longer time periods, errors due to random sources may be reduced and remaining errors can be more representative of biases in the retrieval. The reductions in random errors depend on their correlations in time, however (e.g., Taylor, 1997). When random errors within events are assumed perfectly positively correlated, end-of-event P_ACR accumulations have fractional uncertainties from 1.5 % to 52.4 % (Fig. 12). In actuality, the random error sources likely decorrelate with increasing separation in time. While the scales for these decorrelations are not known, with even a modest amount of decorrelation in the errors the uncertainties are reduced substantially. After applying a negative exponential decorrelation model with a decorrelation scale of 0.5 h to intra-event errors, the fractional uncertainties at the ends of individual events are 1.3 % to 18.8 %. The most significant reductions due to decorrelation occur with the longer-duration events. The end-of-season P_ACR accumulation uncertainties, calculated assuming inter-event uncertainties are uncorrelated, are reduced from 64.9 % for perfectly correlated to 11.8 % for decorrelated intra-event errors.

Figure 12 End-of-event accumulations and uncertainties. The P_ACR accumulation uncertainties are estimated assuming intra-event errors are perfectly correlated (orange) and decorrelated using a negative exponential model with a decorrelation scale of 0.5 h (purple). P_FD12P accumulations (blue-green) are shown for comparison except for those equal to zero, which are omitted. For clarity, the P_ACR accumulations are plotted at ±0.02 h (purple/orange) of their actual durations.

Agreement between observed P_FD12P event accumulations and those from P_ACR generally improves for events with larger accumulations and durations (Fig.12). Of the seven events with accumulations larger than 0.2 mm and durations of 1 h and longer, the P_FD12P accumulations for six fall within or near the uncertainty bounds of the P_ACR accumulations with perfectly correlated errors, while four out of seven are within or near the much narrower bounds for errors with decorrelations. This result is also true for the season as a whole. For the duration of 26.3 h and accumulation of 5.05 mm from P_ACR, the difference compared to the P_FD12P seasonal accumulation of 5.77 mm is −12.6 %. The difference is similar to the P_ACR accumulation uncertainty of 11.7 % for decorrelated errors.

4.2 Information content

The optimal estimation results allow easy calculation of a number of metrics that quantify retrieval performance in terms of information content (Rodgers, 2000; Shannon and Weaver, 1949). These include the averaging kernel matrix

$$\begin{array}{}\text{(22)}& \mathbf{A}={\left({\widehat{\mathbf{K}}}^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\widehat{\mathbf{K}}+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\widehat{\mathbf{K}}}^{\mathsf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\widehat{\mathbf{K}},\end{array}$$

the Shannon information content

$$\begin{array}{}\text{(23)}& H={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\log }_{\mathrm{2}}\left|{\mathbf{S}}_a{\widehat{\mathbf{S}}}_x^{-\mathrm{1}}\right|,\end{array}$$

and the degrees of freedom for signal

$$\begin{array}{}\text{(24)}& d_S=\mathrm{Tr}\left(\mathbf{A}\right).\end{array}$$

Briefly, the diagonal values of A indicate the degree to which the corresponding retrieved state variables are determined by the observations (values nearer 1) versus by the a priori (values nearer 0). H measures how well the observations serve to narrow the possible retrieved states in comparison to the a priori state. Its value can be interpreted as describing the binary bits of resolution of the observing system (L'Ecuyer et al., 2006). d_S quantifies the number of independent quantities that are determined by the observations. See Rodgers (2000) for a more complete discussion in the context of retrieval theory.

For the ACR retrievals, values for H vary between 0.4 and 1.2 (Fig. 13), indicating that the measurements resolve between 1.3 and 2.3 distinct states. Values for d_s show that the retrieval produces somewhat less than one independent piece of information that is significant compared to the measurement and forward model uncertainties. Figure 13c, d show the diagonal elements of A. While the element relevant to λ, A[log (λ)], is consistently positive, the element for N_o, A[log (N₀)], is near zero and is at times negative. These results show that log (λ) is moderately to strongly constrained by the reflectivity observation, while log (N₀) is largely dependent on the a priori constraint.

Figure 13 Distributions of information content metrics for the ACR retrieval.

The size distribution plays a significant role in determining the values of these metrics. Information content H increases as the distribution narrows (Fig. 14a). The increase in H accompanies a substantial increase in the magnitude of the sensitivity of the forward model to log (λ) (Fig. 14b). In contrast, the sensitivity to log (N₀) has a constant value of 10 owing to the reflectivity in dBZ_e being a linear function of log (N₀) (and so is not shown in Fig. 14). This increased sensitivity to log (λ) allows the observed reflectivity to better constrain the retrieved state, particularly the value of log (λ). As a result, A[log (λ)] increases from 0.4 to 0.95 as λ increases (Fig. 14c). The behavior of A[log (N₀)] (Fig. 14d) is quite different. The values are small and are positive for small values of λ but become negative as λ increases. This behavior results from the positive a priori correlation between log (λ) and log (N₀) and the opposing signs of the sensitivities of dBZ_e to these two variables. While the forward model is strongly sensitive to log (λ), its sensitivity to log (N₀) is 3–4 times smaller in magnitude. Consequently, the retrieved value of log (λ) is influenced more strongly by the observations, while the retrieved value of log (N₀) is influenced more by the a priori estimate of the state. This difference is reflected in panels (c) and (d) of Fig. 14.

Figure 14 Information content metrics and the forward model Jacobian as functions of λ.

5 Discussion and conclusions

While millimeter-wavelength, single-frequency radar reflectivity observations alone would seem to have limited utility for retrieving snowfall properties, the results herein demonstrate capabilities for quantifying snowfall rate, accumulation, and aspects of the snow PSD. The results were obtained by applying the radar observations to constrain a priori information appropriate to a broad range of snowfall regimes. The results indicate that the approach would provide useful information when applied to observations such as those from satellite-borne radars, which observe a range of snowfall regimes and for which radar observables are limited to reflectivity.

The results demonstrate the ability of the retrieval to produce reliable estimates of snow accumulation, particularly over timescales involving multiple events and more than several hours of snowfall duration, in spite of large uncertainties in retrieved instantaneous snowfall rates. For the C3VP season, the retrieval reproduced the observed accumulation within 13 % at the end of the season. These results were achieved by omitting two particular time periods during which the retrieval's particle property assumptions were likely very inconsistent with the observed snowfall. Without this adjustment, the end-of-season absolute difference was 18.9 %, illustrating the need for adequate discrimination of the precipitation phase in the retrieval process. Keeping in mind that certain a priori assumptions of the retrieval were also sourced from the C3VP observations, these results are probably best viewed as indicating proper function of the retrieval. The time series of seasonal accumulation shows that while the initial fractional differences reach almost 80 %, the differences diminish with time and increasing accumulation, reaching values of less than 5 % within 5 h. These results are partly due to offsetting errors between events; however, for individual events, best agreement between the observed and retrieved snow accumulations were achieved for events that were longer in duration and produced more substantial accumulations. The observed accumulations for these events were mostly near or within the tighter uncertainty bounds produced by a decorrelating error model applied to the retrieved accumulation. The modest decorrelation used in the model produces uncertainties in the retrieved event accumulations of only 1 % to 20 %. Thus, despite large instantaneous snowfall rate uncertainties for these single-frequency, millimeter-wavelength retrievals, retrieved rates can be expected to prove of value for quantifying accumulations over events, months, seasons, and longer.

Uncertainties in instantaneous retrieval-estimated snowfall rate are dominated by uncertainties in the retrieved state (the uncertainties in the estimated PSD), followed by uncertainties in particle model parameters and, to a lesser extent, the uncertainties in the fall-speed model. Uncertainties due to the assumption of an exponential PSD form are negligible. There is a degree of ambiguity here. The uncertainties in the particle model parameters contribute to the uncertainties in the estimated snowfall rate due to the appearance of the mass term in Eq. (18) but also contribute to uncertainties in the retrieved state. We treat these as independent contributions to the uncertainty. There is likely some covariance that could reduce total uncertainties, but this is not addressed in the treatment of snowfall rate uncertainty presented here.

Retrieval performance, quantified in terms of information content metrics, is determined by the sensitivity of the observations to the desired state vector, the uncertainties assessed for the forward models and measurements, and the explicit assumptions about the uncertainty in the a priori knowledge of the state. For W-band modeled reflectivities in dB, the magnitudes of sensitivities for log (λ) are 3–4 times those for log (N₀), and sensitivities are opposite in sign. This contributes to log (λ) being better constrained by the retrieval than is log (N₀). The consequences of these sensitivities are described more fully in Appendix A. To the extent that process information can be gleaned from changes in the slope parameter over time or space, the retrieval may be useful for process analyses when more direct observations of PSD are not available.

Model-measurement uncertainties are dominated by uncertainties in the particle model parameters (e.g., the coefficients and exponents of the mass- and area-dimension relationships, Table 2), and it is the uncertainties in mass parameters that are the most substantial contributor (Wood et al., 2015). For these near-surface observations, contributions to uncertainties in W-band radar reflectivity from shape, the assumption of an exponential form for the PSD, and the discrete-truncated form of the integrations over size distribution were not significant. For longer wavelength radars that might be used in similar applications (e.g., the MRR or KAZR), shape uncertainty will likely be even smaller due to less prevalent non-Rayleigh scattering.

Table 2 Contributions to uncertainties in forward-modeled and observed reflectivity.

These baseline results suggest several avenues for improving such single-frequency, radar reflectivity-based snowfall retrievals. Improved constraints on snow PSD parameters, through either reduced a priori uncertainties or better observational constraints, are paramount. For ground- or aircraft-based observations, ancillary measurements of snow PSDs can improve the a priori constraints. For retrievals from satellite-borne radar where such measurements are not available, the a priori state is given by more broadly applicable relationships for PSD parameters like those presented here. To the extent that a priori states for specific snowfall regimes might have smaller uncertainties, knowledge of regime-specific PDFs for snow PSD parameters would improve retrieval results provided the correct regime can be diagnosed by the retrieval. Coincident dual-frequency radar observations may also provide improved constraints on the snow PSD parameters (Liao et al., 2005; Matrosov, 2011) but among current satellite-borne instruments, the CPR is single-frequency, and while the GPM DPR provides dual-frequency observations, the DPR sensitivities limit observations to heavier snowfall (Skofronick-Jackson et al., 2019) and implementation of dual-frequency snowfall retrieval has proven difficult (Iguchi et al., 2018). Finally, model-measurement uncertainties can be reduced by reducing uncertainties in particle mass estimates. This may require a more synergistic approach in which improved PSD information is coupled with additional observations such as Doppler velocity to better constrain the assumed particle model used in the retrieval, e.g., moving toward the approach used by Wood et al. (2014, 2015) with ground-based observations. The methods presented here, easily adaptable to other observing systems providing multiple frequency or collocated Doppler velocity observations, provide the basis from which such improvements can be tested and evaluated.

Appendix A: Retrieval interpretation

To interpret the behavior of the retrieval, we refer to the discussion of the information content metrics (Sect. 4.2). The small values for A[log (N₀)] indicate its value is determined primarily by the a priori information and the negative signs do not fit the normal paradigm used to explain the A matrix. Their explanation reveals details of the significant behavior of this retrieval. In the application of the retrieval to a single radar bin, the value of A[log (N₀)] is given by

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\mathbf{A}\left[\log \left(N_{\mathrm{0}}\right)\right]& =\left[s^{\mathrm{2}}\left(\log \left({\widehat{N}}_{\mathrm{0}}\right)\right){\left({\displaystyle \frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(N_{\mathrm{0}}\right)}}\right)}^{\mathrm{2}}+\right.\\ & \left.s\left(\log \left({\widehat{N}}_{\mathrm{0}}\right),\log \left(\widehat{\mathit{\lambda }}\right)\right)\left({\displaystyle \frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(N_{\mathrm{0}}\right)}}\right)\left({\displaystyle \frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(\mathit{\lambda }\right)}}\right)\phantom{{\left({\displaystyle \frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(N_{\mathrm{0}}\right)}}\right)}^{\mathrm{2}}}\right]\\ & {\left[s_y^{\mathrm{2}}\left({\mathrm{dBZ}}_e\right)\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where the carets indicate retrieved values. In the first set of brackets on the right side, the sign of the first term is clearly positive, while that of the second term depends on the signs of the covariance and the two partial derivatives, which are the elements of the Jacobian of the forward model. As was shown earlier (Fig. 14), $\frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(N_{\mathrm{0}}\right)}$ is positive while $\frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(\mathit{\lambda }\right)}$ is negative. The covariance for the retrieved state changes very little from the a priori covariance, which is positive and represents a substantial correlation between log (λ) and log (N₀). This second term, then, is negative and as the magnitude of $\frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(\mathit{\lambda }\right)}$ increases, the sign of A[log (N₀)] changes from positive to negative.

These terms represent competing influences on the retrieved value of log (N₀). These competing influences arise from the a priori covariance and from the Jacobian of the forward model. The positive covariance requires that a positive adjustment in log (λ) be accompanied by a positive adjustment in log (N₀). In contrast, the Jacobian terms have differing signs. If the difference between the observed and forward model reflectivity calls for a positive adjustment to log (λ), the corresponding adjustment to log (N₀) would be negative.

Figure A1 shows this process schematically. The size distribution that represents the initial state is shown by the soid line. Assuming that the forward modeled reflectivity for this state overestimates the observed reflectivity (a positive error), two responses are possible: log (λ) could be increased, narrowing the distribution; and log (N₀) could be decreased, reducing the amplitude of the distribution. Absent the covariance between log (λ) and log (N₀), the retrieval would apply both adjustments, likely giving more weight to the adjustment of log (λ) because of the stronger sensitivity of the forward model to that variable. These adjustments are represented by the heavy arrows labeled δlog (N₀) and δlog (λ). Because of the positive covariance between log (N₀) and log (λ), however, an increase in log (λ) produces an opposing response that increases log (N₀), shown by the upward-pointing heavy arrow. The resulting size distribution is shown by the dashed line.

For small λ (broad distributions), the magnitude of $\frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(\mathit{\lambda }\right)}$ is relatively small, so the covariance-driven adjustment is small and does not overcome the initial reduction in log (N₀). In these cases, log (N₀) decreases in response to a positive error in the modeled reflectivity. This net response is consistent with the sensitivity of the forward model to log (N₀) and A[log (N₀)] is positive. For large λ (narrower distributions), the magnitude of $\frac{\partial {\mathrm{dBZ}}_e}{\partial \log \left(\mathit{\lambda }\right)}$ is larger. The covariance-driven adjustment is larger also and does overcome the initial reduction in log (N₀). As a result, log (N₀) increases in response to the positive error in the modeled reflectivity. Since this net response opposes the sensitivity of the forward model, A[log (N₀)] is negative.

The combination of the strong positive covariance between log (N₀) and log (λ) and the comparatively weak sensitivity of the reflectivity to log (N₀) limits the behavior of the retrieval. For narrower distributions, the retrieval is prevented from simultaneously increasing log (λ) and decreasing log (N₀) in response to a positive error in reflectivity. The opposing behavior, decreasing log (λ) and increasing log (N₀) in response to a negative error in reflectivity, is also restricted. While correct in a climatological sense since log (λ) and log (N₀) are positively correlated, in nature there are likely scenes for which such responses would give a more accurate retrieval. This reasoning demonstrates how other measurements, specifically those with better sensitivity to log (N₀), would benefit the retrieval.

Figure A1 Schematic illustration of the retrieval process. The solid line represents the initial state of the retrieval, while the dashed line shows the adjusted state assuming the initial state overestimates the observed reflectivity. The arrows labeled δlog (λ) and δlog (N₀) show the expected responses of the retrieval based on the sensitivities of the forward model. The arrow labeled s(log (N_o),log (λ)) shows the response due to positive covariance between λ and N₀.

Appendix B: Particle model

The properties here are for the particle shape denoted as “B8pr-30” from Wood et al. (2015), an idealized eight-arm branched spatial particle. Values for the parameters of the mass- and area-dimension power functions are

$$\begin{array}{l}\ln \left(\mathit{\alpha }\right)=-\mathrm{5.723},\\ \mathit{\beta }=\mathrm{2.248},\\ \ln \left(\mathit{\gamma }\right)=-\mathrm{1.379},\\ \mathit{\sigma }=\mathrm{1.813},\end{array}$$

with error covariance matrix

$${\mathbf{S}}_b=\begin{array}{cccc}\mathrm{0.592}& \mathrm{0.212}& \mathrm{0.090}& \mathrm{0.023}\\ \mathrm{0.212}& \mathrm{0.142}& \mathrm{0.011}& \mathrm{0.007}\\ \mathrm{0.090}& \mathrm{0.011}& \mathrm{0.335}& \mathrm{0.103}\\ \mathrm{0.023}& \mathrm{0.007}& \mathrm{0.103}& \mathrm{0.046}\end{array}.$$

These values are appropriate for use with particle size D_M in centimeters, mass in grams and area in square centimeters. The radar backscatter and extinction cross sections are given in Table B1 versus particle size.

Table B1 Backscatter and extinction properties for the snow particle model.