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 radar-only 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}{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}$. The attenuated reflectivity can be written as

where A_dB is the two-way 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 lower-frequency signals are generally attenuated less. In the case of snowfall, the W-band signal can be significantly attenuated, the Ka band much less so, and the Ku band is practically unattenuated by the snowflakes. Thus, the Ku-band radar reflectivity can be used to correct the Ka- and W-band 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 Ku-band 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 root-mean-square (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 ITU-R P.676-11 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 two-way 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 mixed-phase 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₀ 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 large-sized 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² or m is that in the Rayleigh scattering regime (D≪λ) σ_bsc and σ_sca are proportional to m², and σ_abs is proportional to m (Bohren and Huffman, 1983). It follows that the normalized cross sections are roughly constant at the small-particle 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₀, Λ, α 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₀ and Λ, which allows us to compute the integral in Eq. (1) with fixed-point 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₀, Λ 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 dual-wavelength ratios (DWRs), such as simultaneous measurements of Ka–W-band and Ku–Ka-band DWRs, contain information about the size and density of the snowflakes. Following this concept, we form the measurement vector with the Ku-band reflectivity and the Ka–W-band and Ku–Ka-band DWRs. The measurement vector is then

The choice of the Ku-band 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 dual-frequency radar, y consists of the reflectivity from the lowest-frequency radar and the DWR. For single-frequency 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(y|x) 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[x|y]. Multilinear interpolation is used to estimate E[x|y] 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[x|y] 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{‾}}{x}}_i-\mathrm{3}{\mathit{\sigma }}_i$,${\stackrel{\mathrm{‾}}{x}}_i+\mathrm{3}{\mathit{\sigma }}_i$) along each variable, where ${\stackrel{\mathrm{‾}}{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_x|y, 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[x|y] 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, Λ⁻¹ 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 mass-weighted mean diameter

Similarly, the total number concentration of snowflakes may be a more meaningful quantity than N₀. 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.

5.1 3 December 2015

Figure 1The paths of the flights used in Sect. 5.1 (a) and 5.2 (b). The darker sections of the paths show the flight data used in this study (the rest of the measurements were discarded for the lack of useful data). The time stamps (UTC) denote the beginning and end of each flight and the beginning and end of the data that were used. The gray background shows the outline of the Olympic Peninsula, with Vancouver Island to the north.

Download

The first of the two cases that we examined together with NPOL data took place on 3 December 2015. The APR-3 flight leg started at 16:17:23 UTC over the Olympic Mountains, from where the DC-8 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.

Figure 2Data from the 3 December 2015 case described in Sect. 5.1. (a) The mass-weighted mean diameter D_m (Eq. 14). (b) The bulk density ρ_bulk (Eq. 16). (c) The ice water content (Eq. 12). (d) The NPOL hydrometeor identification. (e) A scatter plot of D_m and ρ_bulk from (a) and (b), with the red and orange points identifying the data inside the boxes of corresponding colors shown in those panels. (f) The radar reflectivity observed by NPOL.

Download

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 triple-frequency 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 APR-3 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}{}^{\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 three-frequency 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 DC-8 followed a flight path similar to in the previous case (Fig. 1b); APR-3 data collection for the dataset shown here started at 14:53:21 UTC. We collocated two NPOL RHIs to APR-3 coordinates, one pointing toward land and the other toward the ocean. For the ocean-pointing scan, we selected an RHI that is offset by 4^∘ from the optimal collocation with APR-3 in order to better capture a convective plume that was observed by APR-3 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.

Figure 3As Fig. 2, except for the 4 December 2015 case described in Sect. 5.2. The location of NPOL on the flight track is marked by the arrow in panels (a)–(d) and (f).

Download

The large convective plume found by APR-3 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 APR-3 and NPOL appears to be 1–2 km; this is consistent with the 13 m s⁻¹ 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 graupel-containing 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 APR-3. 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.

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 dual-frequency retrievals and two DWRs for the triple-frequency 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 triple-frequency retrieval are similar to those of the dual-frequency retrievals. However, the multifrequency retrievals clearly outperform single-frequency retrievals. The triple- and dual-frequency scatter plots are visually similar for all two- and three-frequency combinations for Λ, and to a lesser extent N₀. The dual-frequency retrieval using the Ka–W bands seems to be limited in its ability to determine the size of large particles (small Λ), presumably because the dual-frequency ratio saturates at large sizes, while the Ku–Ka-band retrieval suffers from a similar problem with small particles. The Ku–W-band retrieval and the triple-frequency retrieval do not suffer from this problem. Meanwhile, the single-frequency retrievals all have poor sensitivity to N₀. Ku- and Ka-band single-frequency retrievals have some sensitivity to Λ for small particles, while the W-band retrieval also cannot discern this parameter particularly well. None of the retrievals perform adequately with α, although the multifrequency retrievals, especially the triple-frequency retrieval, permit considerably more variation in the values of that parameter: α is almost constant with the single-frequency retrievals, while its relative standard deviation is about 60 % in the triple-frequency 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 single-frequency 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 three-frequency retrieval seems to yield a modest but fairly consistent improvement over the dual-frequency results. These, like with the in situ data comparison, are clearly better than the single-frequency results for all parameters, although the differences for α, W_ice and ρ_bulk are less pronounced.

Figure 5The average posterior retrieval errors of the logarithms of microphysical variables with different combinations of radar frequencies. The data from the 4 December 2015 case (Sect. 5.2) are used in this figure.

Download

The errors in the single-frequency retrievals are all similar; the W band seems to have somewhat smaller errors for W_ice and N₀ while the Ku band is slightly better with the particle size. Notably, the a posteriori errors for the single-frequency retrievals are not much smaller than the a priori errors of Std_a[ln N₀]=2.45, Std_a[ln Λ]=0.83 and Std_a[ln α]=1.13, which emphasizes the poor information content in the single-frequency retrievals. Regardless, with W_ice the single-frequency retrievals perform nearly as well as the multifrequency ones, consistent with what was shown in the comparison to in situ values. None of the dual-frequency options are significantly better than the others, either, although the Ku–Ka-band configuration underperforms the Ka–W-band and Ku–W-band configurations in retrievals of N₀ and N_T, and to a lesser extent W_ice. The Ka–W- and Ku–W-band 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 triple-frequency plots of D_m and especially ρ_bulk compared to the dual-frequency plots. The Ka–W band dual-frequency retrieval appears to capture the plume found by the triple-frequency approach, albeit with a more subdued signal; the other two dual-frequency configurations miss the plume altogether. Consistent with the results of other comparisons shown in this section, the dual-frequency plots capture more detail than the single-frequency plots. This is especially striking for the plots of ρ_bulk, in which the single-frequency 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 single-frequency 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.

Figure 6The root-mean-square changes in the microphysical parameters in response to changes in the prior. The change in the prior is indicated on the left side of each row. The data are from the 4 December 2015 case (Sect. 5.2).

Download

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₀. 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₀ or lnα induce considerably larger changes in the results. Thus, the triple-frequency 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 single-frequency 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 β.

Figure 7The root-mean-square changes in the microphysical parameters in response to changes in the mass–dimensional exponent β. The standard assumption of this paper, β=2.1, is used as the baseline. The value of β is indicated on the left side of each row. The analysis is based on the 4 December 2015 case (Sect. 5.2).

Download

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₀ 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.