the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Retrieval of snowflake microphysical properties from multifrequency radar observations
Jussi Leinonen
Matthew D. Lebsock
Simone Tanelli
Ousmane O. Sy
Brenda Dolan
Randy J. Chase
Joseph A. Finlon
Annakaisa von Lerber
Dmitri Moisseev
We have developed an algorithm that retrieves the size, number concentration and density of falling snow from multifrequency radar observations. This work builds on previous studies that have indicated that threefrequency radars can provide information on snow density, potentially improving the accuracy of snow parameter estimates. The algorithm is based on a Bayesian framework, using lookup tables mapping the measurement space to the state space, which allows fast and robust retrieval. In the forward model, we calculate the radar reflectivities using recently published snow scattering databases. We demonstrate the algorithm using multifrequency airborne radar observations from the OLYMPEX–RADEX field campaign, comparing the retrieval results to hydrometeor identification using groundbased polarimetric radar and also to collocated in situ observations made using another aircraft. Using these data, we examine how the availability of multiple frequencies affects the retrieval accuracy, and we test the sensitivity of the algorithm to the prior assumptions. The results suggest that multifrequency radars are substantially better than singlefrequency radars at retrieving snow microphysical properties. Meanwhile, triplefrequency radars can retrieve wider ranges of snow density than dualfrequency radars and better locate regions of highdensity snow such as graupel, although these benefits are relatively modest compared to the difference in retrieval performance between dual and singlefrequency radars. We also examine the sensitivity of the retrieval results to the fixed a priori assumptions in the algorithm, showing that the multifrequency method can reliably retrieve snowflake size, while the retrieved number concentration and density are affected significantly by the assumptions.
 Article
(2400 KB) 
Supplement
(1293 KB)  BibTeX
 EndNote
Atmospheric ice formation and growth processes have a major impact on the Earth's radiative balance and on the hydrological cycle. Ice clouds and snowfall occur nearly everywhere, as ice processes occur at high altitudes even in areas where freezing temperatures at the surface are rare (Field and Heymsfield, 2015; Mülmenstädt et al., 2015). Ice clouds have also long been a challenge for weather and climate models (Waliser et al., 2009). Improving the microphysics schemes, which describe nucleation of small ice crystals and their transformation into precipitationsized particles, is also currently an active area of model development in which conceptually new schemes have been recently introduced (Harrington et al., 2013; Morrison and Milbrandt, 2015).
Observational data are needed to evaluate the representation of ice and snow in models. While direct measurements of ice particle properties can be made in situ, such measurements only produce limited samples and are difficult and expensive to make, especially when surface observations are not possible and aircraftbased measurements are needed. Remotesensing instruments are able to sample far larger volumes. Radars, in particular, can make rangeresolved measurements and thus map the vertical structure of the ice cloud–precipitation column. However, the interpretation of radar signatures of ice particles is subject to uncertainties because the microwave scattering properties of icy hydrometeors depend on their size, shape and structure. These are extremely variable, as deposition growth alone results in diverse and often complicated shapes, and further growth through aggregation and riming adds to the complexity (Lamb and Verlinde, 2011; Pruppacher and Klett, 1997).
Multifrequency radars have emerged as a potential tool for ice microphysics investigations. It has been recognized for a while that snowflake size can be constrained with collocated measurements at two different frequencies (Hogan et al., 2000; Liao et al., 2005; Matrosov, 1993, 1998). More recently, several studies have shown, using detailed numerical scattering simulations and empirical evidence, that triplefrequency measurements provide information on both the size and density of icy hydrometeors (Gergely et al., 2017; Kneifel et al., 2011, 2015; Kulie et al., 2014; Leinonen and Moisseev, 2015; Leinonen and Szyrmer, 2015; Leinonen et al., 2012a; Stein et al., 2015; Yin et al., 2017). The availability of this information has been expected to enable more accurate quantitative estimation of ice water content (IWC) and snowfall rate and to provide a method to remotely distinguish and characterize icy hydrometeor growth processes.
Studies on the triplefrequency signatures of snow have, so far, been mostly limited to numerical and theoretical investigations, as well as empirical studies that demonstrated the plausibility of the concept. Only very recently have databases of snow scattering properties covering a wide range of snow growth processes (Kuo et al., 2016; Leinonen and Szyrmer, 2015; Lu et al., 2016) become available, enabling the development of a versatile radar forward model that can produce the radar signatures of many types of snowflakes. This, together with the expanded availability of collocated triplefrequency measurement datasets from field campaigns, has provided the prerequisites for the development of a practical snowfall retrieval algorithm for triplefrequency radars.
In this paper, we introduce a method for retrieving certain microphysical properties of snow – namely, the number concentration, size and density – from multifrequency radar observations. The algorithm is based on a Bayesian framework and uses radar cross sections from detailed snowflake models that cover a wide range of sizes and densities. In Sect. 2, we describe the algorithm formulation. Section 3 describes the datasets used for demonstrating and evaluating the algorithm, and Sect. 4 describes how the a priori distributions used in the retrieval were derived. In Sect. 5, we investigate case studies of airborne radar data from the Olympic Mountain Experiment–ACE Radar Definition Experiment 2015 (OLYMPEX–RADEX'15) coordinated by NASA and compare the retrieval results to groundbased polarimetric radar observations. Section 6 describes comparisons to collocated in situ measurements. Section 7 presents statistical analyses of the sensitivity of the algorithm to the number of frequencies available and to the a priori assumptions. Finally, we discuss the implications of the results and summarize the study in Sect. 8.
2.1 Physical basis
The objective of a radar retrieval algorithm for snowfall is to provide the best estimate of the microphysical properties of the snowflakes based on the received radar signals. The unattenuated equivalent radar reflectivity factor Z_{e} for a given wavelength λ is
where σ_{bsc}(D) is the backscattering cross section as a function of the maximum diameter D, N(D) is the particle size distribution and K_{w} is the dielectric factor defined as ${K}_{w}=({n}_{\mathrm{w}}^{\mathrm{2}}\mathrm{1})/({n}_{\mathrm{w}}^{\mathrm{2}}+\mathrm{2})$, where n_{w} is the complex refractive index of liquid water assumed at a reference temperature and frequency.
The attenuation of the radar signal must be accounted for in radaronly retrieval algorithms. The attenuated reflectivity at distance r from the radar is given by
where σ_{ext} is the extinction cross section. The resulting reflectivity is usually expressed in logarithmic units of decibels relative to Z (dBZ), defined by
where ${Z}_{\mathrm{0}}=\mathrm{1}\phantom{\rule{0.125em}{0ex}}{\mathrm{mm}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}$. The attenuated reflectivity can be written as
where A_{dB} is the twoway specific attenuation, that is, the attenuation in decibels per unit length.
It was shown as early as Hitschfeld and Bordan (1954) that weather radar attenuation correction is subject to mathematical instabilities that can lead to small errors multiplying in a positive feedback loop. Namely, overestimation of attenuation in one radar range bin leads to overcompensation in all subsequent bins away from the radar, causing overestimation of the precipitation signal, which in turn leads to further overestimation of the attenuation. In multifrequency radars, the lowerfrequency signals are generally attenuated less. In the case of snowfall, the Wband signal can be significantly attenuated, the Ka band much less so, and the Ku band is practically unattenuated by the snowflakes. Thus, the Kuband radar reflectivity can be used to correct the Ka and Wband signals in a stable manner.
We use a technique similar to Kulie et al. (2014) for attenuation correction: we draw samples from the a priori distribution (described in Sect. 4), calculate both the Kuband reflectivity and the specific attenuation at the Ka or W band for each sample, and fit a function between the reflectivity and the attenuation. We found that a linear function between Z_{dB} and ln A_{dB} fits the relationship well. We validated this approach by computing attenuation afterwards from the retrieved microphysical values; the rootmeansquare (RMS) difference in the total attenuation, calculated over all bins in the case shown in Sect. 5.2, is only 0.27 dB, so this approximate approach to attenuation correction seems to work adequately.
Attenuation also results from atmospheric gases and from supercooled liquid water. The gaseous attenuation was calculated and corrected for with the ITUR P.67611 model (ITU, 2016), using radio sounding data for the temperature, pressure and humidity required by the model. The gaseous attenuation varies spatially since it is dependent on water vapor, but the error introduced by this is likely small given that the maximum twoway gaseous attenuation in the cases analyzed in this study is only 1.1 dB at 94 GHz (W band), and much less at the lower frequencies. However, supercooled liquid water found in mixedphase clouds can cause significant radar attenuation. However, the radar echo of the supercooled water is very weak because of the small size of the drops, making it practically impossible, using radar signals alone, to detect supercooled water coexisting with ice. Thus, we do not correct for attenuation caused by supercooled water, while acknowledging its role as a potential error source.
In order to manage the dimensionality of the problem, the microphysical properties of the snowflakes must be parameterized. We utilize two common assumptions for this. First, we assume that the particle size distribution (PSD) follows the exponential distribution
where N_{0} and Λ are the intercept and slope parameters, respectively. Although gamma distributions, and other forms that introduce additional parameters, are sometimes used, the exponential distribution has been found to describe snowflake size distributions well (Heymsfield et al., 2008; Sekhon and Srivastava, 1970). We also found it to be a good match to the in situ airborne size distribution measurements used in this study (see Sect. 6). Therefore, we find it preferable over more complicated alternatives. Second, we assume that the mass of snowflakes is given as a function of the diameter as
as has been commonly done in microphysics literature (Pruppacher and Klett, 1997). In the following section, we explain how these assumptions are used to compute the radar reflectivities.
2.2 Forward model
The forward model in an inversion algorithm is responsible for calculating the measurements that correspond to a given state vector – in our case, the radar reflectivity at each wavelength given the microphysical parameters. The simulation of radar reflectivity from snowflakes whose diameters are comparable to or larger than the wavelength is known to require calculations that account for the internal structure of the snowflake (Botta et al., 2011; Petty and Huang, 2010; Tyynelä et al., 2011). Recently, such calculations have been used for a wide variety of model snowflakes in order to establish databases of scattering properties. We chose a combination of two such datasets as the basis of our forward model: the rimed snowflakes of Leinonen and Szyrmer (2015) and the OpenSSP database of Kuo et al. (2016). The dataset of Leinonen and Szyrmer (2015) covers a wide range of snowflake densities, but due to the relatively coarse resolution of the volume elements, it mostly contains moderate and largesized snowflakes. The Kuo et al. (2016) dataset was used to augment the set of snowflakes used by the forward model at small sizes, D<1 mm.
While there have been considerable recent advances on the problem of modeling snowflakes produced by different ice processes and calculating their scattering properties, the abundance of available snowflake models leads to another question: which set of snowflakes should be used by the forward model in a particular situation? We use an approach that does not force us to select any one dataset. Instead, the scattering properties are given as a function of mass and size: σ(D,m), where σ can be one of σ_{bsc}, σ_{sca} or σ_{abs}, the last two being the scattering and absorption cross sections, respectively, with ${\mathit{\sigma}}_{\mathrm{ext}}={\mathit{\sigma}}_{\mathrm{sca}}+{\mathit{\sigma}}_{\mathrm{abs}}$. The function σ(D,m) is constructed by organizing all model snowflakes from the combined scattering database into bins by D and m; we use 128×128 logarithmically spaced bins to cover the range of diameters and masses found in the dataset. For each bin, we compute the average of ${\mathit{\sigma}}_{\mathrm{norm}}\equiv \mathit{\sigma}/{m}^{\mathit{\gamma}}$, where γ=2 for the backscattering and scattering cross sections, and γ=1 for the absorption cross section. The reason for the normalization by m^{2} or m is that in the Rayleigh scattering regime (D≪λ) σ_{bsc} and σ_{sca} are proportional to m^{2}, and σ_{abs} is proportional to m (Bohren and Huffman, 1983). It follows that the normalized cross sections are roughly constant at the smallparticle limit. To smoothen the binned values, the samples used in the averaging are weighted using a Gaussian function of the distance from the bin center, with a standard deviation of 0.15 for both lnD and lnm. A continuous function of the form
is then formed by interpolation among the bin centers. Not all bins have snowflakes in them; for those we are unable to do the averaging and instead assign the scattering properties to zero. This means that the limits of the coverage of the snowflake database in the (D,m) space are effectively assumed to be the limits of the natural variability in snowflakes. While this is not exactly true, the combined database does cover a wide range of microphysical processes. The assumption that the cross section goes to zero (as opposed to, for instance, extrapolating it) outside the coverage area also effectively truncates the integrals in Eqs. (1) and (2).
With a method to calculate the cross sections as a function of D and m, it is relatively straightforward to compute radar reflectivities from the microphysical inputs. As can be seen from the previous section, the input parameters for the forward model are N_{0}, Λ, α and β. We start with a fixed set of 1024 logarithmically spaced integration points that span the interval [D_{min},D_{max}]. The parameters α and β are used to find the corresponding masses using Eq. (6). The cross section for each integration point is then found from the lookup table using interpolation. The cross sections are multiplied with the size distribution determined by N_{0} and Λ, which allows us to compute the integral in Eq. (1) with fixedpoint numerical integration.
2.3 Retrieval
A radar retrieval algorithm needs to invert Eqs. (1) and (2) such that an input of ${Z}_{\mathrm{e}}^{\prime}$ at one or more wavelengths yields the properties of N(D) and m(D). The inversion is unavoidably inexact, as the wide variety of snowflake number concentrations, size distributions and densities leads to a variability too large to constrain with a few radar reflectivities. The retrieval must be performed in a probabilistic sense, deriving the most likely solution from the possible alternatives, using the prior information about snowflake properties as a constraint.
The retrieval problem is commonly stated as finding a state vector x that explains a given measurement vector y. The formulation of the state vector depends on which variables are chosen for retrieval and which ones are simply assumed. In our experimentation with different combinations, we found that the most stable solution was to retrieve N_{0}, Λ and α. The β parameter was fixed at 2.1. While β varies in nature, many experimental and modeling studies (Delanoë et al., 2014; Erfani and Mitchell, 2017; Leinonen and Moisseev, 2015; Mascio and Mace, 2017; Mascio et al., 2017; Mitchell et al., 1990; Moisseev et al., 2017; Pruppacher and Klett, 1997; Westbrook et al., 2004) have found exponents near this value for various types of snowflakes; we will examine the sensitivity of the results to this assumption in Sect. 7.3. We retrieve the logarithm of each microphysical parameter because the dynamic ranges of the retrieved values are large and because using the logarithmic values makes the forward model more linear; this was examined analytically for the simpler case of cloud water retrieval by Leinonen et al. (2016). The state vector then becomes
In our multifrequency radar retrieval algorithm, the most straightforward way to formulate the measurement vector would be to use each of the three radar reflectivities. However, earlier studies (Kneifel et al., 2011; Leinonen and Szyrmer, 2015) have shown that combinations of dualwavelength ratios (DWRs), such as simultaneous measurements of Ka–Wband and Ku–Kaband DWRs, contain information about the size and density of the snowflakes. Following this concept, we form the measurement vector with the Kuband reflectivity and the Ka–Wband and Ku–Kaband DWRs. The measurement vector is then
The choice of the Kuband reflectivity is somewhat arbitrary, as any of the three bands could be used, but the Ku band does benefit from that band being the least attenuated of the three. In studies in which we omit one of the radar bands, instead operating with a dualfrequency radar, y consists of the reflectivity from the lowestfrequency radar and the DWR. For singlefrequency retrievals, y simply contains Z_{dB} at the single band.
The measurement vector must be accompanied by an error estimate, which should include not only the radar instrument error but also the error due to the forward model assumptions. In our case, the latter includes the errors due to the assumptions of an exponential size distribution, a fixed mass–dimensional exponent β and the orientation distributions assumed in the scattering databases. The extent of these errors is difficult to quantify, but their effect should be similar on each collocated radar frequency: for example, the radar cross section will increase with increasing particle size for all frequencies, and thus the errors in radar reflectivity at different frequencies will partially cancel out when computing the DWRs. This suggests that the DWRs can be assumed to have smaller errors than the absolute value of the reflectivity. Accordingly, we assign 3 dB of error standard deviation for the absolute value of the radar reflectivity and 1 dB for each of the DWRs.
In atmospheric remote sensing, the inversion problem is often solved using optimal estimation (Rodgers, 2000). This is a Bayesian method that assumes that x and y are jointly distributed according to the multivariate normal distribution and which is solved using optimization methods. We found this technique to be problematic for our retrieval, partly due to the limited and discrete nature of the snowflake scattering database used in the forward model. The optimization in OE often converged to local minima, especially near the extreme values supported by the snowflake database, introducing sudden changes to the retrieved values.
Despite the shortcomings of OE, a Bayesian approach was still desirable in order to constrain the retrieved microphysical parameters. We found that the retrieval can be performed in a robust way through a global calculation of the expected value of the state x given a measurement y. This is given by
where p(y) is the marginal probability of y, p(yx) is the conditional probability of a measurement y given a state x and p(x) is the a priori probability of x, described in detail in Sect. 4. This approach is slightly different from the common strategy of finding the most likely solution given the prior and the measurement: That method aims to find the mode of the conditional distribution; ours determines the mean.
Using Eq. (10), we can construct a lookup table that maps discrete values of y to the corresponding expected values E[xy]. Multilinear interpolation is used to estimate E[xy] for values of y that fall between the discrete values used in the table. The errors associated with the discretization can be reduced to be arbitrarily small by making the intervals between the values finer. In the studies presented here, the lookup table for E[xy] ranged between 0 and 35 dBZ for Z_{dB,Ku}, between −2 and 14 dB for DWR_{Ka∕W}, and between −2 and 9 dB for DWR_{Ku∕Ka}, with 0.25 dB discretization for each dimension. The integral in Eq. (10) was computed by evaluating the integrand at approximately 10 000 discrete points, which were distributed uniformly across a finite search space spanning (${\stackrel{\mathrm{\u203e}}{x}}_{i}\mathrm{3}{\mathit{\sigma}}_{i}$,${\stackrel{\mathrm{\u203e}}{x}}_{i}+\mathrm{3}{\mathit{\sigma}}_{i}$) along each variable, where ${\stackrel{\mathrm{\u203e}}{x}}_{i}$ is the prior mean of the ith variable in x, and σ_{i} is its prior standard deviation. Making the discretization finer than this did not seem to change the retrieval results significantly in our case, although we encourage those using this approach for other problems to establish the appropriate discretization for their problem.
Error estimates for the retrieved values can be computed using the same technique. The error covariance matrix of the state given an observation, S_{xy}, can be computed as
where “⊗” is the outer product. $\mathbf{E}[\mathit{x}\otimes \mathit{x}\mathit{y}]$ can be evaluated using a lookup table and interpolation in the same manner as explained for E[xy] above.
The method described above allows the state and its covariance to be retrieved robustly and very quickly, with only a table lookup and an interpolation needed for each measurement. This comes at the cost of a relatively expensive initialization of the tables before the retrieval is started. However, with our parameters for the discretization, this only took about 1 min on a modern laptop computer with no parallelization, so it does not present a major computational burden.
2.4 Derived variables
The results of the retrieval are the parameters of Eq. (8), but for further analysis of the results, it is useful to derive other variables that are important for microphysics or more intuitively understood by end users. Perhaps most importantly, the IWC (denoted by W_{ice}), which expresses the ice mass in a unit volume of air, is given by
Consistently with the calculation of the scattering properties, we set m(D)=0 in the integral (and other integrals in this section) where no snowflake samples are available for the (D,m) combination. If this truncation is not used, the assumptions of Eqs. (5) and (6) give W_{ice} in the simple form
where Γ is the gamma function.
When discussing the snowflake size, Λ^{−1} gives the average diameter for the untruncated exponential size distribution, but it is often clearer and more convenient to state the diameter that contributes most to the IWC. This is the massweighted mean diameter
Similarly, the total number concentration of snowflakes may be a more meaningful quantity than N_{0}. This is given simply by
Also, the density of the snowflakes depends on the diameter, but a bulk density for the snowflake ensemble can be computed by dividing the IWC by the volume spanned by the enclosing spheres of the snowflakes in a unit volume:
We use this definition for simplicity; a somewhat higher density would be obtained by using the volume of the enclosing spheroid or ellipsoid in the integral in the denominator, but the shape of this ellipsoid is in general dependent on D and m, which would complicate the calculation.
The quantities in Eqs. (12)–(16) are nonlinear functions of the state x, and consequently estimating their errors is not completely straightforward. Since our algorithm returns a probability distribution function (PDF) for x, we can obtain statistically valid error estimates by computing the standard deviation of a quantity over the PDF. This can be estimated quickly with Gauss–Hermite quadratures; see Appendix A.
The main source of data that we use to demonstrate the triplefrequency retrieval is from the Airborne Third Generation Precipitation Radar (Sadowy et al., 2003) flown onboard the NASA DC8 aircraft during the OLYMPEX–RADEX experiment, which took place around the Olympic Mountains of Washington State, USA, in late 2015 (Houze et al., 2017). The RADEX involvement in this field campaign was intended specifically to assess the the capabilities of multifrequency radar observations for satellite remote sensing of precipitation processes. APR3 acquired simultaneous measurements at three frequencies: 13.4 GHz (Ku band), 35.6 GHz (Ka band) and 94.9 GHz (W band). APR3 is a scanning polarimetric cloudprofiling radar with Doppler capability. With a vertical resolution of 30 m, it provides highresolution 3D measurements of clouds and precipitation. OLYMPEX was the first time it was deployed in its triplefrequency configuration.
We investigated the ability of the triplefrequency algorithm to identify snowfall processes qualitatively by comparing the results to collocated groundbased dualpolarization radar observations. These observations were made by the NASA SBand DualPolarimetric Radar (NPOL), which was deployed 2 km from the coast at 47.277^{∘} N, 124.211^{∘} W, 157 m above mean sea level (m.s.l.). The NPOL scanning strategy interleaved planned position indicator scans (PPIs) with a series of highresolution rangeheight indicator (RHI) sector scans to the west over the ocean and to the east over the Quinault River valley (Houze et al., 2017). During OLYMPEX, the NASA DC8 aircraft frequently flew directly along NPOL RHI azimuths, making it relatively straightforward to collocate with the nadirpointing scans from APR3. We collocated NPOL data to the APR3 radar coordinates using the Python ARM Radar Toolkit (Helmus and Collis, 2016) by first identifying RHI scans whose time and direction coincided with the APR3 overpass, then copying data from the nearest NPOL bin to each APR3 bin. We used two variables from NPOL: the radar reflectivity and the hydrometeor identification (HID) product (Dolan and Rutledge, 2009). The latter uses fuzzy logic to assign the most likely hydrometeor class to each radar bin based on temperature and the radar reflectivity and polarimetric parameters. We use this product to provide independent estimates of the type of icy hydrometeors and compare them to the microphysical properties retrieved by our algorithm.
During the OLYMPEX campaign, the University of North Dakota Citation aircraft often flew in the same area as the NASA DC8. Typically, the Citation flew at lower altitudes than the DC8, and consequently there are many data points where the Citation measurements are collocated with the APR3. A total of 16 cases from OLYMPEX were analyzed. The APR3 gate closest to the Citation is found using a kdimensionaltree search algorithm. The Citation measured the PSD using the 2DS (Stereo) Probe (Lawson et al., 2006) in the range of $\mathrm{225}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}\le D<\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{mm}$ and the HighVolume Particle Spectrometer ($\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{mm}\le D\le \mathrm{3.25}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}$). To eliminate shattered artifacts created from ice crystals colliding with the probe housing, antishattering tips are used in conjunction with the University of Illinois Oklahoma Optical Array Probe Processing Software (Jackson et al., 2014). In addition to the optical array probes, the Citation also carried a Nevzorov probe (Korolev et al., 1998) to measure bulk total water content.
The groundbased observations of snowfall microphysics used to derive the a priori distribution were gathered at the Hyytiälä Forestry Field Station (61.845^{∘} N, 24.287^{∘} E, 150 m above mean sea level) of the University of Helsinki, Finland, during the Biogenic Aerosols – Effects on Clouds and Climate (BAECC) campaign (Petäjä et al., 2016) and the following winter of 2014–2015. The weather conditions during BAECC and the following winter were mostly mild, and most of the snowfall observations were collected at temperatures above $\mathrm{4}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$. Both aggregation and riming occurred frequently during the measurement period (Moisseev et al., 2017). The PSDs were measured with a video disdrometer, the Snowflake Video Imager (Newman et al., 2009), as a function of the diskequivalent diameter (the diameter of a disk with the projected area of the particle image). The mean PSD was calculated for every 5 min period. The resolution of the SVI is 0.1 mm, although in practice, the smallest diskequivalent diameter used in the computations was approximately 0.2 mm. The PSD was divided into 120 bins with a bin size of 0.2 mm; the highest bin is for diameters larger than 26.0 mm. A linear scaling factor between the diskequivalent diameter and the maximum diameter was determined by analyzing SVI images of snowflakes from each case and utilized to give the PSD as a function of maximum diameter. The mass retrievals were obtained by combining SVI observations with a collocated precipitation gauge. Based on the particle fall velocity and shape measurements provided by the SVI, the masses of individual falling snow particles were estimated with hydrodynamic theory (Mitchell and Heymsfield, 2005; von Lerber et al., 2017). The mass–dimensional relation in the form of Eq. (6) was determined for every 5 min with mass as a function of maximum diameter and with a linear regression fit in the log scale.
We also used balloon sounding data to support the analysis of the case studies. These data were derived from publicly available operational soundings launched daily at 00:00 and 12:00 UTC from Quillayute, Washington, near the area where the radar measurements took place.
Bayesian retrievals depend on the availability of a priori data. We based our a priori values on two sources of in situ data: the Citation dataset from OLYMPEX and the groundbased measurements from BAECC. Both of these datasets can be used to derive the N_{0}, Λ and α parameters. For both datasets, N_{0} and Λ can be derived from the binned PSDs. The α parameter can be derived by fitting a curve defined by Eq. (6) to the mass as a function of diameter; this is included in the BAECC data, in which the mass was derived from the snowflake fall velocity (von Lerber et al., 2017). In calculating α from the BAECC dataset, we fixed β to 2.1, consistent with the assumptions in the retrieval algorithm. For the Citation data, mass is not directly available as a function of diameter, but W_{ice} is estimated with the Nevzorov probe and thus α can be roughly estimated using Eq. (13).
For the purposes of demonstrating the algorithm, we based the a priori distribution used in this study on a combination of the two datasets, taking an equal number of samples from each for a total N≈6000. We recognize that this is an imperfect solution, and a further analysis using these and other datasets should be conducted to establish a priori distributions suitable for remotesensing retrievals of snowfall under various atmospheric conditions. Doing this rigorously will likely require an entire study of its own.
The analysis resulted in means of $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}{N}_{\mathrm{0}}}=\mathrm{15.4}$, $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}\mathrm{\Lambda}}=\mathrm{7.50}$, and $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}\mathit{\alpha}}=\mathrm{2.30}$ and standard deviations of Std[ln N_{0}]=1.67, Std[ln Λ]=0.52 and Std[ln α]=0.69. Because the two datasets cannot be expected to cover the entire natural distribution of these parameters, basing the a priori distribution on them would likely result in an overly restrictive prior. To compensate for this, we increase the standard deviations given above by a factor of 1.5, acknowledging that this choice is somewhat arbitrary. The correlation matrix of x derived from the datasets is
from which the a priori covariance matrix can be computed as
where D is a diagonal matrix with the standard deviations of x on the diagonal. The resulting distribution, used as the prior in all retrievals in this study, is then given by the mean x_{a} and covariance S_{a}:
In Sect. 7.2 we examine the sensitivity of the results to the choice of prior.
We assume that the a priori distribution is multivariate normal. Given the limited scope of the datasets used to derive the prior distribution in this study, we cannot rigorously test this assumption, but the choice is motivated by probabilistic arguments that the normal distribution is the most natural choice for an unknown distribution (Jaynes, 2003). Global distributions for microphysical quantities have also often been found to be lognormal (Kedem and Chiu, 1987; Leinonen et al., 2012b), meaning that the distributions of their logarithms (we use the logarithmic values in the state vector) are normal. Thus a multivariate normal distribution is a reasonable assumption for this study, although larger datasets should be analyzed in this manner in order to derive appropriate global priors.
5.1 3 December 2015
The first of the two cases that we examined together with NPOL data took place on 3 December 2015. The APR3 flight leg started at 16:17:23 UTC over the Olympic Mountains, from where the DC8 flew towards the coast, passing directly over the NPOL site. A map of the flight path is shown in Fig. 1a. The case consisted primarily of prefrontal stratiform precipitation; see Houze et al. (2015a) for details. We only used data from regions above the melting layer, which we identified just below 3 km in altitude based on the the radar bright band; this also agrees with the 0 ^{∘}C isotherm of 2.85 km in the 12:00 UTC balloon sounding from nearby Quillayute, Washington.
The retrievals from the case are shown in Fig. 2a–e. On the left side of Fig. 2a–c, an orange box delineates a column in which D_{m} increases significantly with decreasing altitude, accompanied by a rapid decrease in ρ_{bulk}. Together, these changes point to the onset of aggregation, which results in rapid growth of snowflakes accompanied by a decrease in density as single ice crystals stick together to form aggregates, whose density decreases as a function of size. The transition can also be seen in Fig. 2e, in which orange dots denote the data points from the orange box in Fig. 2a–c.
The transition from ice crystals to aggregates is also detected at around 5 km in altitude, 20–45 km on the distance scale, by both the triplefrequency retrieval, which shows a sudden increase in D_{m} (Fig. 2a), and by NPOL, which identifies a change in the hydrometeor type at roughly the same altitude. According to NPOL HID, the hydrometeors above this altitude consist mostly of a mixture of ice crystals and aggregates, while the hydrometeors below it are identified as aggregates. While the altitude at which aggregation initiates appears to be similar between NPOL and our retrieval, small discrepancies are to be expected because the APR3 observations are not perfectly simultaneous with the NPOL scan. The time difference ranges from 4 min at the beginning of the observations shown in Fig. 2 to 14 min at the end. Further evidence for aggregation is provided by sounding data, which indicate a temperature between −15 and $\mathrm{12}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$ in the layer at 5.0–5.5 km in altitude, a common temperature range for the onset of aggregation driven by dendritic growth of ice crystals at these temperatures (Bailey and Hallett, 2009; Lamb and Verlinde, 2011).
Another interesting feature found in this case is denoted by the red boxes in Fig. 2a–c. In this region, the retrieved microphysical variables indicate moderately sized snowflakes with relatively high ρ, which suggests that rimed snowflakes occur in the area. The data points located within this box are shown in red in the scatter plot of Fig. 2e, which confirms these attributes. It is interesting to note that the red region and the bottom of the orange region have similar IWCs, but the sizes and densities are very different. NPOL also detects some graupel in this region, which suggests that the threefrequency retrieval detects snowflake riming and graupel formation. In the following case, we further explore this capability.
5.2 4 December 2015
On 4 December 2015, precipitation originated mostly from postfrontal convection following the passage of the front on the previous day (Houze et al., 2015b). The DC8 followed a flight path similar to in the previous case (Fig. 1b); APR3 data collection for the dataset shown here started at 14:53:21 UTC. We collocated two NPOL RHIs to APR3 coordinates, one pointing toward land and the other toward the ocean. For the oceanpointing scan, we selected an RHI that is offset by 4^{∘} from the optimal collocation with APR3 in order to better capture a convective plume that was observed by APR3 but had moved before being scanned by NPOL 2 min later. This shifted the location of the scan by only 500 m at the distance of the plume. The sounding data and the radar bright band both placed the melting level at around 1.3 km, lower than on the previous day. We again only used data points located above the melting layer.
The large convective plume found by APR3 in this case is marked with a red box in Fig. 3a–c. As with Fig. 2, the data points from this box are denoted with red dots in Fig. 3e. In this case, the data points from the plume are particularly distinct from the rest of the joint distribution of D_{m} and ρ_{bulk}, indicating moderately large particles with high density, characteristic of graupel. NPOL also indicates a similarly sized plume of graupel in this region. The time separation of the scans in the region to the right of NPOL in Fig. 3 is only 2 min, so it seems likely that the same plume was captured by both radars. The spatial shift between the plumes observed by APR3 and NPOL appears to be 1–2 km; this is consistent with the 13 m s^{−1} wind speed measured by the sounding at 3 km in altitude, which translates to a 1.5 km distance over 2 min.
On the left side of NPOL, another graupelcontaining region is denoted by an orange box. This region is also accompanied by an NPOL detection of graupel in the vicinity. The time separation in this region was longer, between 4 and 8 min, so the plume had more time to move away from the vertical cross section before being observed by APR3. Regardless, the two radars agree on location of the plume to within 2 km and on its height to within 0.5 km.
Our retrieval and the NPOL HID also seem to be in reasonably good agreement regarding the transition from ice crystals to aggregates. Both indicate the presence of ice crystals (i.e., small, relatively dense hydrometeors) at higher altitudes and aggregates at lower altitudes (below approximately 4 km), with the transition point varying considerably within this case. Both products also identify the presence of smaller particles at 2–3 km in altitude in the region, around 20 km on the horizontal scale.
As described in Sect. 3, the UND Citation aircraft gathered particle probe measurements simultaneously with the NASA DC8 radar observations during OLYMPEX. This resulted in a set of collocated radar and in situ data. The retrieval algorithm was run using the collocated and attenuationcorrected radar reflectivity values. The retrieved microphysical quantities were then compared to those measured in situ. The availability of variables from the in situ dataset is somewhat limited: while the number and projected sizes of the ice particles can be measured quite accurately using the imaging probes, the twodimensional nature of the imager limits the accuracy of the maximum dimension as this must be estimated from a projection of the particle. The snowflake masses are also difficult to determine. The bulk IWC can be estimated with the Nevzorov probe, but its inlet is only 8 mm in diameter, which causes it to underestimate IWC when the maximum particle size exceeds approximately 4 mm (Korolev et al., 2013). Unfortunately, the cases with large snowflakes are where one would expect the largest benefits from multifrequency methods because of the stronger resonance effects involved in scattering. Thus, this limitation of the Nevzorov probe somewhat diminishes its value in validating the retrievals. While the Citation measurements do not give the masses of individual particles, α can be estimated from Eq. (13) if the IWC given by the Nevzorov probe is assumed to be correct.
To filter out outliers and poor collocations, we applied two filters. First, to ensure an acceptably accurate collocation between the two measurements, the time separation between them was required to be less than 2 min. Second, for adequate sampling, the total number concentration N_{T} was required to be more than 10^{3} m^{−3}. These criteria successfully removed most outliers that we found in the unfiltered comparisons.
The comparisons of the retrievals against the in situ values are shown on the top row of Fig. 4 (the same analysis run with reduced frequencies, shown on the other rows of Fig. 4, is discussed in Sect. 7.1). The figures show that retrievals of the slope parameter Λ compare considerably better to the in situ values than do the retrievals of the intercept parameter N_{0}, which in turn are better than those of the mass–dimensional factor α. The Λ parameters agree well throughout the range of values (for ln Λ, rootmeansquare error is RMSE[ln Λ]=0.41, bias is bias[ln Λ]=0.023 and correlation is Cor[ln Λ]=0.70), showing that particle sizing can be carried out reliably using the multifrequency retrieval. N_{0} is also quite well matched (Cor[ln N_{0}]=0.56), but the relative errors are much larger than for Λ (RMSE[ln N_{0}]=3.01, $\mathrm{bias}\left[\mathrm{ln}\phantom{\rule{0.125em}{0ex}}{N}_{\mathrm{0}}\right]=\mathrm{0.73}$). The α parameters are poorly matched between the two datasets, although the retrieval produces some variation in this parameter. In any case, one should be skeptical of the α comparison as the in situ values have been derived from the Nevzorov probe data, which suffers from the abovementioned problems, and using Eq. (13), which is an approximation. Furthermore, fixing the β parameter may further exacerbate the problem with estimating α.
The retrieved IWCs W_{ice} correspond quite well to the in situ values (RMSE[ln W_{ice}]=0.72, bias[ln W_{ice}]=0.30, Cor[ln W_{ice}]=0.67). Interestingly, the IWC, which is a function of N_{0} and α, appears to be better retrieved than either of those parameters. Opposite errors in N_{0} and α, seen in their respective scatter plots, suggest that their retrieval errors compensate for each other in a way that allows W_{ice} to be constrained better than either N_{0} or α alone. This is also supported by the correlation matrix of the retrieval errors, for which the error correlation between ln N_{0} and ln α is −0.30 on average. The red dots in Fig. 4 correspond to larger ice particles, where the Nevzorov probe might be prone to underestimation. However, there does not appear to be a significant difference in W_{ice} between the small and large particles. However, the large particles stand out in the α scatter plot, where they are clearly the worst match between the in situ and retrieved values.
7.1 Sensitivity to the number of frequencies
In the assessment of a multifrequency algorithm, one interesting question is what are the benefits of introducing additional frequencies? To evaluate this, we reran the analysis of Sect. 6 with subsets of the frequencies used in the full analysis. We examined all the possible combinations of available bands, always using the lowest frequency for the absolute reflectivity, combined with the DWRs that were available (one DWR for dualfrequency retrievals and two DWRs for the triplefrequency retrieval).
The scatter plots of the in situ and retrieved microphysical parameters are shown in Fig. 4. These plots suggest that the results of the triplefrequency retrieval are similar to those of the dualfrequency retrievals. However, the multifrequency retrievals clearly outperform singlefrequency retrievals. The triple and dualfrequency scatter plots are visually similar for all two and threefrequency combinations for Λ, and to a lesser extent N_{0}. The dualfrequency retrieval using the Ka–W bands seems to be limited in its ability to determine the size of large particles (small Λ), presumably because the dualfrequency ratio saturates at large sizes, while the Ku–Kaband retrieval suffers from a similar problem with small particles. The Ku–Wband retrieval and the triplefrequency retrieval do not suffer from this problem. Meanwhile, the singlefrequency retrievals all have poor sensitivity to N_{0}. Ku and Kaband singlefrequency retrievals have some sensitivity to Λ for small particles, while the Wband retrieval also cannot discern this parameter particularly well. None of the retrievals perform adequately with α, although the multifrequency retrievals, especially the triplefrequency retrieval, permit considerably more variation in the values of that parameter: α is almost constant with the singlefrequency retrievals, while its relative standard deviation is about 60 % in the triplefrequency results, indicating that the retrieval algorithm is confident enough in the signal to estimate α as something other than the a priori mean. The results for α should be interpreted skeptically because of the issues with the derivation of α, as explained in Sect. 6. The singlefrequency retrievals appear to constrain W_{ice} much better than they constrain any of the individual microphysical parameters.
Another way to evaluate the sensitivity to the number of frequencies is to examine the a posteriori errors reported by the algorithm itself. These errors, derived from the 4 December 2015 case, are shown in Fig. 5 for the different frequency combinations. According to the error estimate from the algorithm, the threefrequency retrieval seems to yield a modest but fairly consistent improvement over the dualfrequency results. These, like with the in situ data comparison, are clearly better than the singlefrequency results for all parameters, although the differences for α, W_{ice} and ρ_{bulk} are less pronounced.
The errors in the singlefrequency retrievals are all similar; the W band seems to have somewhat smaller errors for W_{ice} and N_{0} while the Ku band is slightly better with the particle size. Notably, the a posteriori errors for the singlefrequency retrievals are not much smaller than the a priori errors of Std_{a}[ln N_{0}]=2.45, Std_{a}[ln Λ]=0.83 and Std_{a}[ln α]=1.13, which emphasizes the poor information content in the singlefrequency retrievals. Regardless, with W_{ice} the singlefrequency retrievals perform nearly as well as the multifrequency ones, consistent with what was shown in the comparison to in situ values. None of the dualfrequency options are significantly better than the others, either, although the Ku–Kaband configuration underperforms the Ka–Wband and Ku–Wband configurations in retrievals of N_{0} and N_{T}, and to a lesser extent W_{ice}. The Ka–W and Ku–Wband configurations are nearly equally good.
We have additionally created plots of the microphysical parameters shown in Fig. 3 using each of the frequency combinations found in Fig. 5. Due to the large number of plots resulting from this analysis, these plots are not shown here, but can be found in Figs. S1–S21 of the Supplement accompanying this article. A notable feature of these plots is the higher level of detail and wider range of variation found in the triplefrequency plots of D_{m} and especially ρ_{bulk} compared to the dualfrequency plots. The Ka–W band dualfrequency retrieval appears to capture the plume found by the triplefrequency approach, albeit with a more subdued signal; the other two dualfrequency configurations miss the plume altogether. Consistent with the results of other comparisons shown in this section, the dualfrequency plots capture more detail than the singlefrequency plots. This is especially striking for the plots of ρ_{bulk}, in which the singlefrequency retrievals appear to always give nearly the same density. In contrast to D_{m} and ρ_{bulk}, W_{ice} has only small differences, and similar levels of detail between the singlefrequency and multifrequency retrievals. This is again similar to the findings in Fig. 4.
7.2 Sensitivity to prior assumptions
In order to examine the sensitivity of the results of the retrieval algorithm to the prior assumptions, we ran the case of 4 December 2015 with shifted prior means. We changed the mean of each variable in the state vector x, one at a time, by ±1 standard deviation of that variable. The results are shown in Fig. 6. The results are consistent with the retrievals in the sense that a shift in the prior of a variable causes a smaller shift of the same sign in the a posteriori value of that variable.
The effects on other variables from adjusting the prior of one variable are not straightforward to interpret. These are connected in a complicated way due to the significant a priori correlations among the different variables, as well as the necessity of explaining the observed reflectivities with other parameters when one of them is shifted. The dependencies are clearly not linear. The shifts in the prior also interact with the limits of the scattering database, which further complicates the interpretation. The IWC is the most sensitive to the prior of ln N_{0}. The results are the least sensitive to the prior assumption of ln Λ, indicating that ln Λ is very well constrained by the observations. Changes to the priors of either ln N_{0} or lnα induce considerably larger changes in the results. Thus, the triplefrequency algorithm is clearly still somewhat dependent on the a priori assumptions, although the changes in the posterior values are much smaller than the corresponding changes in the prior, showing that the radar signal constrains them quite effectively.
In Figs. S22–S28, we repeat this analysis with the reduced frequencies. These clearly show the increasing dependence on the prior assumptions with fewer available frequencies. Again, the difference between triple and dual frequency is fairly modest, while the singlefrequency retrievals shift much more in response to changes in the prior.
7.3 Sensitivity to mass–dimensional exponent
The most significant fixed parameter in the retrieval is the exponent β of the mass–dimensional relationship (Eq. 6). Similar to Sect. 7.2, we carried out an analysis of the sensitivity of the retrieval results to the choice of β. We used the value usually adopted in this paper, β=2.1, as the reference and compared the results obtained with β=1.9, β=2.3 and β=2.5 to the reference retrieval. The values were chosen based on exponents found in the literature for single crystals, aggregate snowflakes and rimed particles (Mitchell et al., 1990); higher exponents such as those close to 3.0 often found for graupel (Heymsfield and Kajikawa, 1987; Locatelli and Hobbs, 1974) were not tested because the distribution of particles in the scattering databases does not support such high exponents well. The results are shown in Fig. 7. This figure is similar to Fig. 6, but we have omitted the changes in the mass–dimensional prefactor α because this parameter does not have a physical meaning independent of β.
The changes in the retrieval results for different values of β exhibit patterns similar to those resulting from the change in prior values: The parameters corresponding to number concentration (N_{0} and N_{T}) and density (ρ_{bulk}) are the most sensitive to the assumptions. Meanwhile, parameters related to particle size (Λ and D_{m}) and, to a lesser extent, the IWC W_{ice} are less affected by changes in β. The changes in retrieved parameters with changing β can be substantial, suggesting that a good estimate of β is important for quantitatively correct retrievals. However, the changes are predictable and reasonable, which suggests that the algorithm is robust and can function with different values of β without major problems. A notable exception to the predictable behavior is that of W_{ice}, whose retrieved value increases in response to both increase and decrease in β from 2.1.
In this study, we described and evaluated an algorithm for snow microphysical retrievals using multifrequency radar measurements. The probabilistic method is based on direct application of Bayes' theorem using lookup tables. We examined the capabilities and limitations of the retrieval algorithm using data from the OLYMPEX–RADEX measurement campaign, comparing the results to groundbased radar measurements from the NASA NPOL radar and to in situ measurements from the UND Citation aircraft, both of which were collocated with the APR3 measurements. We also examined the sensitivity of the algorithm to various assumptions used in its formulation.
The results indicate that, at least for the retrieval approach presented here, triplefrequency radar retrievals provide modest benefits over dualfrequency retrievals of snowfall properties. The probabilistic error estimates from the triplefrequency retrievals are generally only slightly smaller than those from dualfrequency retrievals, but closer examination of the retrieved values shows that the triplefrequency approach produces more detailed retrievals with higher degrees of variability than the dualfrequency retrievals. The triplefrequency method can also determine particle size throughout the range of snowflake sizes studied here, avoiding problems with some of the dualfrequency methods with sizing either small or large particles. Multifrequency retrievals significantly outperform those using only one frequency, and none of the three dualfrequency configurations studied (Ka–W, Ku–Ka and Ku–Wbands) appear to be decisively better than the others, although the Ka–W band combination was found to have more sensitivity to the snowflake density than the Ku–Ka or Ku–Wband combinations. Similarly, we found the relative performances of Ku, Ka and Wband singlefrequency retrievals to be approximately equal. Thus, information content analysis appears to suggest that multifrequency radars are preferable to singlefrequency radars in snowfall retrievals, but it does not provide much insight into the exact choice of frequencies; this choice should probably be more dependent on other factors such as achievable sensitivity and resolution, the importance of attenuation, and cost.
The triplefrequency technique appears to be useful at identifying graupel, that is, ice particles that are heavily rimed and thus considerably denser than most aggregate snowflakes, providing a sufficient signal for the triplefrequency retrieval to detect. This was confirmed in this study with the comparison to polarimetric observations with the NPOL groundbased radar. Globally, graupel occurs in relatively rare events that represent only a small fraction of snow cases, and consequently graupel events do not impact the statistics much. However, graupel (and hail, which is even denser) can have a substantial societal impact where it occurs, and thus detecting it can be valuable even though it only occurs in a small percentage of icy precipitation. Detecting graupel plumes, together with accurate snowflake size determination elsewhere in a precipitating region, can also shed light on the processes involved in the formation of graupel. These plumes are usually small in their horizontal extent, of the order of 1 km, requiring a fairly high spatial resolution in the radars used to detect them, which can be challenging to achieve if multifrequency radars are considered for satellite applications.
Despite the improvements in retrieval precision in multifrequency retrievals, the retrieved results are still dependent on the assumptions regarding the a priori distribution of the retrieved microphysical parameters, as well as the mass–dimensional exponent β. Different retrieved parameters have widely different sensitivities to the assumptions: the retrieved snow particle size changes only modestly in response to changes to the prior and to β, indicating that the size can be retrieved robustly with the multifrequency method. In contrast, the retrieved number concentration and density are much more sensitive to the assumptions and therefore potentially susceptible to retrieval errors caused by inaccurate prior data. Therefore, it is still vital to constrain the algorithm using in situ measurements that provide not only the size and number concentration of snowflakes but also their mass–dimensional scaling parameters α and β. Later versions of the algorithm should include β as a retrievable parameter and incorporate it in the multivariate prior so that the retrieval errors originating from the uncertainty of β can be properly quantified.
The findings of this study concern the retrieval accuracy of multifrequency radars and do not address their other potential benefits. For instance, multifrequency radars can utilize lowerfrequency channels (e.g., Ku band) to penetrate deeper into precipitation, particularly heavy rain that can attenuate higher frequencies (e.g., W band) heavily enough to block detection altogether. Conversely, higherfrequency radars can generally be made more sensitive, allowing detection in regions below the sensitivity thresholds of lowfrequency bands. These benefits should be considered together with the retrieval performance when decisions about instrument specifications are made; see, e.g., Leinonen et al. (2015) for a quantitative assessment of retrieval capabilities of a potential spaceborne triplefrequency radar.
This work builds on earlier experimental and modeling results that suggested that triplefrequency radars can be used to constrain snowflake habits and examines this capability in practice with a prototype retrieval algorithm. Based on the experience gained in this study, we can identify two requirements for future research that need to be fulfilled in order to use such an algorithm in an operational setting. First, the snowflake scattering database, while more extensive than those previously available, is still limited in its scope, and its coverage of snowflake sizes, densities and habits should be expanded in order to support the forward model in all scenarios. Second, the a priori distributions used in the retrievals in this study are based on relatively few data points. An abundance of in situ data from ice clouds and snowfall currently exists as a result of many ground and aircraftbased field campaigns; analyses of the data from these are needed to support retrieval algorithm development by providing representative a priori distributions of snowfall properties. The substantial cross correlations found in this study among the snow microphysical properties (Eq. 17) emphasize the need for a multivariate analysis of these datasets.
The APR3 data files can be downloaded from the OLYMPEX data repository at https://doi.org/10.5067/GPMGV/OLYMPEX/APR3/DATA201 (Durden and Tanelli, 2018), and the NPOL data from https://doi.org/10.5067/GPMGV/OLYMPEX/NPOL/DATA301 (Wolff et al., 2017). The Citation data are available at https://github.com/dopplerchase/Chase_et_al_2018 (last access: 1 October 2018), maintained by Randy J. Chase (email: randyjc2@illinois.edu). The BAECC campaign data are available at https://github.com/dmoisseev/SnowRetrievals20142015 (last access: 1 October 2018). The sounding data can be obtained from the University of Wyoming collection at http://weather.uwyo.edu/upperair/sounding.html (last access: 1 October 2018). The retrieval results, used to generate the plots, are available in numerical form from Jussi Leinonen (email: jussi.s.leinonen@jpl.caltech.edu).
Consider a scalar Q(x) that is a function (not necessarily a linear function) of the vector x of normally distributed random variables, whose probability distribution p(x) is given by the mean 〈x〉 and the covariance S. For example, Q can be ln W_{ice} or the logarithm of any variable introduced in Sect. 2.4. Then, a probabilistic error estimate is given by the standard deviation
where the expectation, denoted by 〈⋅〉, is taken over the PDF of x. The expectation can be estimated efficiently using a Gauss–Hermite quadrature. For a threevariable x (generalization to other numbers of variables is straightforward), the expectation 〈Q〉 is obtained as follows:
where

V is a matrix whose columns contain the normalized eigenvectors of S,

Λ is a diagonal matrix containing the corresponding eigenvalues of S,

x_{GH} and w_{GH} are the points and weights of a Gauss–Hermite quadrature that gives the approximation
$$\begin{array}{}\text{(A5)}& \underset{\mathrm{\infty}}{\overset{\mathrm{\infty}}{\int}}\mathrm{exp}({x}^{\mathrm{2}})f\left(x\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}x\approx \sum _{i=\mathrm{1}}^{N}{w}_{\mathrm{GH},i}\phantom{\rule{0.125em}{0ex}}f\left({x}_{\mathrm{GH},i}\right),\end{array}$$where the approximation is exact if f is a polynomial of at most degree 2N−1; x_{GH} and w_{GH} can be found in many tables (Beyer, 1987) and in scientific software packages (Oliphant, 2007).
〈Q^{2}〉 can also be estimated using the above method, thus giving the error estimate when substituted into Eq. (A1). This is derived by computing the Gauss–Hermite quadrature for the standard multivariate normal distribution with zero mean and identity covariance, then mapping the quadrature points to the corresponding points in the distribution of x.
The supplement related to this article is available online at: https://doi.org/10.5194/amt1154712018supplement.
JL planned the study, formulated and implemented the retrieval algorithm, and performed the analysis presented in this paper. He also led the preparation of this article, with contributions from all authors. MDL advised on the algorithm formulation and data analysis. ST and OOS calibrated and quality controlled the APR3 data and provided support for using them; ST also participated in the collection of the APR3 data during OLYMPEX–RADEX. BD processed the NPOL data, provided advice on their use and participated in the NPOL operations during OLYMPEX. RJC and JAF analyzed and processed the Citation data, collocated them with APR3, and advised on the comparisons to the retrievals. AvL and DM coordinated the collection of the BAECC in situ data and processed them for use in this study.
The authors declare that they have no conflicts of interest.
We thank the two anonymous reviewers for their constructive comments. The
research of Jussi Leinonen, Matthew D. Lebsock, Simone Tanelli and Ousmane O. Sy was carried out at the Jet Propulsion
Laboratory (JPL), California Institute of Technology, under contract with
NASA. The work of Jussi Leinonen and Matthew D. Lebsock was supported by the NASA AerosolCloudEcosystem
and CloudSat missions under RTOP WBS 103930/6.1 and 103428/8.A.1.6,
respectively. Jussi Leinonen was partly funded under subcontract 1559252 from JPL to
UCLA. Simone Tanelli and Ousmane O. Sy acknowledge support from the GPM GV program and the ACE
Science Working Group funding for the acquisition and initial processing of
APR3 data, and support from the Earth Science U.S. Participating Investigator program for
the detailed analysis of Wband Doppler data. Funding for the research of Brenda Dolan,
Randy J. Chase and Joseph A. Finlon was provided by NASA Precipitation Measurement Missions
grants
NNX16AI11G (BD) and NNX16AD80G (Randy J. Chase and Joseph A. Finlon) under Ramesh Kakar. Dmitri Moisseev
acknowledges the funding received through ERAPLANET, transnational project
iCUPE (grant agreement 689443), funded under the EU Horizon 2020 Framework
Programme, and the Academy of Finland (grant nos. 307331 and 305175). The
research work of Annakaisa von Lerber was funded by EU's Horizon 2020 research and
innovation program (ECHORIZON2020PR700099ANYWHERE).
Edited by: Mark Kulie
Reviewed by: two anonymous referees
Bailey, M. P. and Hallett, J.: A Comprehensive Habit Diagram for Atmospheric Ice Crystals: Confirmation from the Laboratory, AIRS II, and Other Field Studies, J. Atmos. Sci., 66, 2888–2899, https://doi.org/10.1175/2009JAS2883.1, 2009. a
Beyer, W. H.: CRC Handbook of Mathematical Sciences, CRC Press, Boca Raton, Florida, USA, 1987. a
Bohren, C. F. and Huffman, D. R.: Absorption and Scattering of Light by Small Particles, John Wiley & Sons, Inc., New York, USA, 1983. a
Botta, G., Aydin, K., Verlinde, J., Avramov, A. E., Ackerman, A. S., Fridlind, A. M., McFarquhar, G. M., and Wolde, M.: Millimeter wave scattering from ice crystals and their aggregates: Comparing cloud model simulations with Xand Kaband radar measurements, J. Geophys. Res., 116, D00T04, https://doi.org/10.1029/2011JD015909, 2011. a
Delanoë, J. M. E., Heymsfield, A. J., Protat, A., Bansemer, A., and Hogan, R. J.: Normalized particle size distribution for remote sensing application, J. Geophys. Res.Atmos., 119, 4204–4227, https://doi.org/10.1002/2013JD020700, 2014. a
Dolan, B. and Rutledge, S. A.: A theorybased hydrometeor identification algorithm for Xband polarimetric radars, J. Atmos. Ocean. Tech., 46, 1196–1213, https://doi.org/10.1175/2009JTECHA1208.1, 2009. a
Durden, S. L. and Tanelli, S.: GPM Ground Validation Airborne Precipitation Radar 3rd Generation (APR3) OLYMPEX V2, Dataset available online from the NASA EOSDIS Global Hydrology Resource Center Distributed Active Archive Center, Huntsville, Alabama, USA, https://doi.org/10.5067/GPMGV/OLYMPEX/APR3/DATA201, 2018.
Erfani, E. and Mitchell, D. L.: Growth of ice particle mass and projected area during riming, Atmos. Chem. Phys., 17, 1241–1257, https://doi.org/10.5194/acp1712412017, 2017. a
Field, P. R. and Heymsfield, A. J.: Importance of snow to global precipitation, Geophys. Res. Lett., 42, 9512–9520, https://doi.org/10.1002/2015GL065497, 2015. a
Gergely, M., Cooper, S. J., and Garrett, T. J.: Using snowflake surfaceareatovolume ratio to model and interpret snowfall triplefrequency radar signatures, Atmos. Chem. Phys., 17, 12011–12030, https://doi.org/10.5194/acp17120112017, 2017. a
Harrington, J. Y., Sulia, K., and Morrison, H.: A Method for Adaptive Habit Prediction in Bulk Microphysical Models. Part I: Theoretical Development, J. Atmos. Sci., 70, 349–364, https://doi.org/10.1175/JASD12040.1, 2013. a
Helmus, J. J. and Collis, S. M.: The Python ARM Radar Toolkit (PyART), a Library for Working with Weather Radar Data in the Python Programming Language, J. Open Res. Software, 4, e25, https://doi.org/10.5334/jors.119, 2016. a
Heymsfield, A. J. and Kajikawa, M.: An Improved Approach to Calculating Terminal Velocities of Platelike Crystals and Graupel, J. Atmos. Sci., 44, 1088–1099, https://doi.org/10.1175/15200469(1987)044<1088:AIATCT>2.0.CO;2, 1987. a
Heymsfield, A. J., Field, P., and Bansemer, A.: Exponential size distributions for snow, J. Atmos. Sci., 65, 4017–4031, https://doi.org/10.1175/2008JAS2583.1, 2008. a
Hitschfeld, W. and Bordan, J.: Errors Inherent in the Radar Measurement of Rainfall at Attenuating Wavelenghts, J. Meteorol., 11, 58–67, https://doi.org/10.1175/15200469(1954)011<0058:EIITRM>2.0.CO;2, 1954. a
Hogan, R. J., Illingworth, A. J., and Sauvageot, H.: Measuring crystal size in cirrus using 35 and 94GHz radars, J. Atmos. Ocean. Tech., 17, 27–37, https://doi.org/10.1175/15200426(2000)017<0027:MCSICU>2.0.CO;2, 2000. a
Houze Jr., R. A., McMurdie, L., Tanelli, S., Mace, J., and Nesbitt, S.: OLYMPEX Science Summary for 3 December 2015, available at: http://olympex.atmos.washington.edu/archive/reports/20151203/20151203Science_summary.html (last access: 2 February 2018), 2015a. a
Houze Jr., R. A., McMurdie, L., Zagrodnik, J., Duffy, G., Durden, S., and Funk, A.: OLYMPEX Science Summary for 4 December 2015, available at: http://olympex.atmos.washington.edu/archive/reports/20151204/20151204Science_summary.html (last access: 2 February 2018), 2015b. a
Houze Jr., R. A., McMurdie, L. A., Petersen, W. A., Schwaller, M. R., Baccus, W., Lundquist, J. D., Mass, C. F., Nijssen, B., Rutledge, S. A., Hudak, D. R., Tanelli, S., Mace, G. G., Poellot, M. R., Lettenmaier, D. P., Zagrodnik, J. P., Rowe, A. K., DeHart, J. C., Madaus, L. E., and Barnes, H. C.: The Olympic Mountains Experiment (OLYMPEX), B. Am. Meteorol. Soc., 98, 2167–2188, https://doi.org/10.1175/BAMSD160182.1, 2017. a, b
ITU: Recommendation ITUR P.67611: Attenuation by atmospheric gases, International Telecommunications Union, 2016. a
Jackson, R. C., McFarquhar, G. M., Stith, J., Beals, M., Shaw, R. A., Jensen, J., Fugal, J., and Korolev, A.: An Assessment of the Impact of Antishattering Tips and Artifact Removal Techniques on Cloud Ice Size Distributions Measured by the 2D Cloud Probe, J. Atmos. Ocean. Tech., 31, 2567–2590, https://doi.org/10.1175/JTECHD1300239.1, 2014. a
Jaynes, E. T.: Probability Theory: The Logic of Science, Cambridge University Press, Cambridge, UK, 2003. a
Kedem, B. and Chiu, L.: On the lognormality of rain rate, P. Natl. Acad. Sci. USA, 84, 901–905, 1987. a
Kneifel, S., Kulie, M. S., and Bennartz, R.: A triple frequency approach to retrieve microphysical snowfall parameters, J. Geophys. Res., 116, D11203, https://doi.org/10.1029/2010JD015430, 2011. a, b
Kneifel, S., von Lerber, A., Tiira, J., Moisseev, D., Kollias, P., and Leinonen, J.: Observed relations between snowfall microphysics and triplefrequency radar measurements, J. Geophys. Res.Atmos., 120, 6034–6055, https://doi.org/10.1002/2015JD023156, 2015. a
Korolev, A., Strapp, J. W., Isaac, G. A., and Emery, E.: Improved Airborne HotWire Measurements of Ice Water Content in Clouds, J. Atmos. Ocean. Tech., 30, 2121–2131, https://doi.org/10.1175/JTECHD1300007.1, 2013. a
Korolev, A. V., Strapp, J. W., Isaac, G. A., and Nevzorov, A. N.: The Nevzorov Airborne HotWire LWCTWC Probe: Principle of Operation and Performance Characteristics, J. Atmos. Ocean. Tech., 15, 1495–1510, https://doi.org/10.1175/15200426(1998)015<1495:TNAHWL>2.0.CO;2, 1998. a
Kulie, M. S., Hiley, M. J., Bennartz, R., Kneifel, S., and Tanelli, S.: Triple frequency radar reflectivity signatures of snow: Observations and comparisons to theoretical ice particle scattering models, J. Appl. Meteorol. Clim., 53, 1080–1098, https://doi.org/10.1175/JAMCD13066.1, 2014. a, b
Kuo, K.S., Olson, W. S., Johnson, B. T., Grecu, M., Tian, L., Clune, T. L., van Aartsen, B. H., Heymsfield, A. J., Liao, L., and Meneghini, R.: The Microwave Radiative Properties of Falling Snow Derived from Nonspherical Ice Particle Models. Part I: An Extensive Database of Simulated Pristine Crystals and Aggregate Particles, and Their Scattering Properties, J. Appl. Meteorol. Clim., 55, 691–708, https://doi.org/10.1175/JAMCD150130.1, 2016. a, b, c
Lamb, D. and Verlinde, J.: Physics and Chemistry of Clouds, Cambridge University Press, Cambridge, UK, 2011. a, b
Lawson, R. P., O'Connor, D., Zmarzly, P., Weaver, K., Baker, B., Mo, Q., and Jonsson, H.: The 2DS (Stereo) Probe: Design and Preliminary Tests of a New Airborne, HighSpeed, HighResolution Particle Imaging Probe, J. Atmos. Ocean. Tech., 23, 1462–1477, https://doi.org/10.1175/JTECH1927.1, 2006. a
Leinonen, J. and Moisseev, D.: What do triplefrequency radar signatures reveal about aggregate snowflakes?, J. Geophys. Res., 120, 229–239, https://doi.org/10.1002/2014JD022072, 2015. a, b
Leinonen, J. and Szyrmer, W.: Radar signatures of snowflake riming: A modeling study, Earth Space Sci., 2, 346–358, https://doi.org/10.1002/2015EA000102, 2015. a, b, c, d, e
Leinonen, J., Kneifel, S., Moisseev, D., Tyynelä, J., Tanelli, S., and Nousiainen, T.: Evidence of Nonspheroidal Behavior in MillimeterWavelength Radar Observations of Snowfall, J. Geophys. Res., 117, D18205, https://doi.org/10.1029/2012JD017680, 2012a. a
Leinonen, J., Moisseev, D., Leskinen, M., and Petersen, W.: A Climatology of Disdrometer Measurements of Rainfall in Finland over Five Years with Implications for Global Radar Observations, J. Appl. Meteorol. Clim., 51, 392–404, https://doi.org/10.1175/JAMCD11056.1, 2012b. a
Leinonen, J., Lebsock, M. D., Tanelli, S., Suzuki, K., Yashiro, H., and Miyamoto, Y.: Performance assessment of a triplefrequency spaceborne cloudprecipitation radar concept using a global cloudresolving model, Atmos. Meas. Tech., 8, 3493–3517, https://doi.org/10.5194/amt834932015, 2015. a
Leinonen, J., Lebsock, M. D., Stephens, G. L., and Suzuki, K.: Improved Retrieval of Cloud Liquid Water from CloudSat and MODIS, J. Appl. Meteorol. Clim., 55, 1831–1844, https://doi.org/10.1175/JAMCD160077.1, 2016. a
Liao, L., Meneghini, R., Iguchi, T., and Detwiler, A.: Use of dualwavelength radar for snow parameter estimates, J. Atmos. Ocean. Tech., 22, 1494–1506, https://doi.org/10.1175/JTECH1808.1, 2005. a
Locatelli, J. D. and Hobbs, P. V.: Fall speeds and masses of solid precipitation particles, J. Geophys. Res., 79, 2185–2197, https://doi.org/10.1029/JC079i015p02185, 1974. a
Lu, Y., Jiang, Z., Aydin, K., Verlinde, J., Clothiaux, E. E., and Botta, G.: A polarimetric scattering database for nonspherical ice particles at microwave wavelengths, Atmos. Meas. Tech., 9, 5119–5134, https://doi.org/10.5194/amt951192016, 2016. a
Mascio, J. and Mace, G. G.: Quantifying uncertainties in radar forward models through a comparison between CloudSat and SPartICus reflectivity factors, J. Geophys. Res.Atmos., 122, 1665–1684, https://doi.org/10.1002/2016JD025183, 2017. a
Mascio, J., Xu, Z., and Mace, G. G.: The MassDimensional Properties of Cirrus Clouds During TC4, J. Geophys. Res.Atmos., 122, 10402–10417, https://doi.org/10.1002/2017JD026787, 2017. a
Matrosov, S. Y.: Possibilities of cirrus particle sizing from dualfrequency radar measurements, J. Geophys. Res., 98, 20675–20683, https://doi.org/10.1029/93JD02335, 1993. a
Matrosov, S. Y.: A dualwavelength radar method to measure snowfall rate, J. Appl. Meteorol., 37, 1510–1521, https://doi.org/10.1175/15200450(1998)037<1510:ADWRMT>2.0.CO;2, 1998. a
Mitchell, D. L. and Heymsfield, A. J.: Refinements in the Treatment of Ice Particle Terminal Velocities, Highlighting Aggregates, J. Atmos. Sci., 62, 1637–1644, https://doi.org/10.1175/JAS3413.1, 2005. a
Mitchell, D. L., Zhang, R., and Pitter, R. L.: MassDimensional Relationships for Ice Particles and the Influence of Riming on Snowfall Rates, J. Appl. Meteorol., 29, 153–163, https://doi.org/10.1175/15200450(1990)029<0153:MDRFIP>2.0.CO;2, 1990. a, b
Moisseev, D., von Lerber, A., and Tiira, J.: Quantifying the effect of riming on snowfall using groundbased observations, J. Geophys. Res.Atmos., 122, 4019–4037, https://doi.org/10.1002/2016JD026272, 2017. a, b
Morrison, H. and Milbrandt, J. A.: Parameterization of Cloud Microphysics Based on the Prediction of Bulk Ice Particle Properties. Part I: Scheme Description and Idealized Tests, J. Atmos. Sci., 72, 287–311, https://doi.org/10.1175/JASD140065.1, 2015. a
Mülmenstädt, J., Sourdeval, O., Delanoë, J., and Quaas, J.: Frequency of occurrence of rain from liquid, mixed, and icephase clouds derived from ATrain satellite retrievals, Geophys. Res. Lett., 42, 6502–6509, https://doi.org/10.1002/2015GL064604, 2015. a
Newman, A. J., Kucera, P. A., and Bliven, L. F.: Presenting the Snowflake Video Imager (SVI), J. Atmos. Ocean. Tech., 26, 167–179, https://doi.org/10.1175/2008JTECHA1148.1, 2009. a
Oliphant, T. E.: Python for Scientific Computing, Comput. Sci. Eng., 9, 10–20, https://doi.org/10.1109/MCSE.2007.58, 2007. a
Petäjä, T., O'Connor, E. J., Moisseev, D., Sinclair, V. A., Manninen, A. J., Väänänen, R., von Lerber, A., Thornton, J. A., Nicoll, K., Petersen, W., Chandrasekar, V., Smith, J. N., Winkler, P. M., Krüger, O., Hakola, H., Timonen, H., Brus, D., Laurila, T., Asmi, E., Riekkola, M.L., Mona, L., Massoli, P., Engelmann, R., Komppula, M., Wang, J., Kuang, C., Bäck, J., Virtanen, A., Levula, J., Ritsche, M., and Hickmon, N.: BAECC: A Field Campaign to Elucidate the Impact of Biogenic Aerosols on Clouds and Climate, B. Am. Meteorol. Soc., 97, 1909–1928, https://doi.org/10.1175/BAMSD1400199.1, 2016. a
Petty, G. W. and Huang, W.: Microwave Backscatter and Extinction by Soft Ice Spheres and Complex Snow Aggregates, J. Atmos. Sci., 67, 769–787, https://doi.org/10.1175/2009JAS3146.1, 2010. a
Pruppacher, H. R. and Klett, J. D.: Microphysics of clouds and precipitation, Springer, Rordrecht, The Netherlands, 2nd edn., 1997. a, b, c
Rodgers, C. D.: Inverse Methods for Atmospheric Sounding – Theory and Practice, World Scientific Publishing, https://doi.org/10.1142/9789812813718, 2000. a
Sadowy, G. A., Berkun, A. C., Chun, W., Im, E., and Durden, S. L.: Development of an advanced airborne precipitation radar, Microwave J., 46, 84, available at: http://www.microwavejournal.com/articles/3577developmentofanadvancedairborneprecipitationradar (last access: 1 October 2018), 2003. a
Sekhon, R. S. and Srivastava, R. C.: Snow size spectra and radar reflectivity, J. Atmos. Sci., 27, 299–307, https://doi.org/10.1175/15200469(1970)027<0299:SSSARR>2.0.CO;2, 1970. a
Stein, T. H. M., Westbrook, C. D., and Nicol, J. C.: Fractal geometry of aggregate snowflakes revealed by triplewavelength radar measurements, Geophys. Res. Lett., 43, 176–183, https://doi.org/10.1002/2014GL062170, 2015. a
Tyynelä, J., Leinonen, J., Moisseev, D., and Nousiainen, T.: Radar backscattering from snowflakes: comparison of fractal, aggregate and softspheroid models, J. Atmos. Ocean. Tech., 28, 1365–1372, https://doi.org/10.1175/JTECHD1100004.1, 2011. a
von Lerber, A., Moisseev, D., Bliven, L. F., Petersen, W., Harri, A.M., and Chandrasekar, V.: Microphysical Properties of Snow and Their Link to Ze'S Relations during BAECC 2014, J. Appl. Meteorol. Clim., 56, 1561–1582, https://doi.org/10.1175/JAMCD160379.1, 2017. a, b
Waliser, D. E., Li, J.L. F., Woods, C. P., Austin, R. T., Bacmeister, J., Chern, J., Del Genio, A., Jiang, J. H., Kuang, Z., Meng, H., Minnis, P., Platnick, S., Rossow, W. B., Stephens, G. L., SunMack, S., Tao, W.K., Tompkins, A. M., Vane, D. G., Walker, C., and Wu, D.: Cloud ice: A climate model challenge with signs and expectations of progress, J. Geophys. Res.Atmos., 114, D00A21, https://doi.org/10.1029/2008JD010015, 2009. a
Westbrook, C. D., Ball, R. C., Field, P. R., and Heymsfield, A. J.: Universality in snowflake aggregation, Geophys. Res. Lett., 31, L15104, https://doi.org/10.1029/2004GL020363, 2004. a
Wolff, D., Marks, D., Petersen, W. A., and Pippitt, J.: GPM Ground Validation NASA SBand Dual Polarimetric (NPOL) Doppler Radar OLYMPEX V2. Dataset available online from the NASA EOSDIS Global Hydrology Resource Center Distributed Active Archive Center, Huntsville, Alabama, USA, https://doi.org/10.5067/GPMGV/OLYMPEX/NPOL/DATA301, 2017.
Yin, M., Liu, G., Honeyager, R., and Turk, F. J.: Observed differences of triplefrequency radar signatures between snowflakes in stratiform and convective clouds, J. Quant. Spectrosc. Ra., 193, 13–20, https://doi.org/10.1016/j.jqsrt.2017.02.017, 2017. a
 Abstract
 Introduction
 Algorithm
 Data
 A priori assumptions
 Case studies and comparison to NPOL
 Comparison to in situ data
 Sensitivity analysis
 Conclusions
 Data availability
 Appendix A: Fast derivation of error estimates for retrieved quantities
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 Algorithm
 Data
 A priori assumptions
 Case studies and comparison to NPOL
 Comparison to in situ data
 Sensitivity analysis
 Conclusions
 Data availability
 Appendix A: Fast derivation of error estimates for retrieved quantities
 Author contributions
 Competing interests
 Acknowledgements
 References
 Supplement