Snow microphysical retrieval from the NASA D3R radar during ICE-POP 2018

A method is developed to use both polarimetric and dual-frequency radar measurements to retrieve microphysical properties of falling snow. It is applied to the Ku- and Ka-band measurements of the NASA dual-polarization, dual-frequency Doppler radar (D3R) obtained during the International Collaborative Experiments for PyeongChang 2018 Olympic and Paralympic winter games (ICE-POP 2018) field campaign and incorporates the Atmospheric Radiative Transfer Simulator (ARTS) microwave single-scattering property database for oriented particles. The retrieval uses optimal estimation to solve for several parameters that describe the particle size distribution (PSD), relative contribution of pristine, aggregate, and rimed ice species, and the orientation distribution along an entire radial simultaneously. Examination of Jacobian matrices and averaging kernels shows that the dual-wavelength ratio (DWR) measurements provide information regarding the characteristic particle size, and to a lesser extent, the rime fraction and shape parameter of the size distribution, whereas the polarimetric measurements provide information regarding the mass fraction of pristine particles and their characteristic size and orientation distribution. Thus, by combining the dual-frequency and polarimetric measurements, some ambiguities can be resolved that should allow a better determination of the PSD and bulk microphysical properties (e.g., snowfall rate) than can be retrieved from single-frequency polarimetric measurements or dual-frequency, single-polarization measurements.

The D3R ICE-POP retrievals were validated using Precipitation Imaging Package (PIP) and Pluvio weighing gauge measurements taken nearby at the May Hills ground site. The PIP measures the snow PSD directly, and its measurements can be used to derived the snowfall rate (volumetric and water equivalent), mean volume-weighted particle size, and effective density, as well as particle aspect ratio and orientation. Four retrieval experiments were performed to evaluate the utility of different measurement combinations: Ku-only, DWR-only, Ku-pol, and All-obs. In terms of correlation, the volumetric snowfall rate (r=0.95) and snow water equivalent rate (r=0.92) were best retrieved by the Ku-pol method, while the DWR-only method had the lowest magnitude bias for these parameters (−31 % and −8 %, respectively). The methods that incorporated DWR also had the best correlation to particle size (r=0.74 and r=0.71 for DWR-only and All-obs, respectively), although none of the methods retrieved density particularly well (r=0.43 for All-obs). The ability of the measurements to retrieve mean aspect ratio was also inconclusive, although the polarimetric methods (Ku-pol and All-obs) had reduced biases and mean absolute error (MAE) relative to the Ku-only and DWR-only methods. The significant biases in particle size and snowfall rate appeared to be related to biases in the measured DWR, emphasizing the need for accurate DWR measurements and frequent calibration in future D3R deployments.

1 Introduction

Estimation of snowfall rates and other properties from weather radar is made difficult by many of the same challenges that exist for rainfall estimation (primarily, the discrepancy between the sixth-moment dependence of radar reflectivity factor Z and the third- to fourth-moment dependence of precipitation rate R), but additional factors further confound radar retrievals of snow. Whereas the shape and scattering properties of a raindrop depend only on its mass and temperature (e.g., Beard et al., 2010; Ekelund et al., 2020), there is tremendous diversity in ice crystals of a given mass, resulting from the infinite complexity of particle trajectories through differing thermodynamic environments resulting in growth by vapor deposition (e.g., Kuroda and Lacmann, 1982; Chen and Lamb, 1994; Fukuta and Takahashi, 1999), aggregation (e.g., Hosler et al., 1957; Hobbs et al., 1974; Connolly et al., 2012), and riming (e.g., Mitchell et al., 1990; Jensen and Harrington, 2015; Moisseev et al., 2017), as well as ablation by sublimation (e.g., Smith et al., 2009) and melting (e.g., Matsuo and Sasyo, 1981; Leinonen and von Lerber, 2018). All of these processes influence the scattering and aerodynamic properties of these ice particles in ways that can influence the interpretation of radar data (e.g., Hall et al., 1984; Vivekanandan et al., 1994; Bechini et al., 2013; Botta et al., 2013; Thompson et al., 2014). Despite these challenges, since the introduction of weather radar, these instruments have been important tools in gathering information about ice-phase precipitation. Early efforts focused on using Z to estimate the intensity of snow precipitation measurements with assumed ice particle size distribution (PSD) forms (e.g., Marshall and Gunn, 1952). These efforts relied on in situ ground measurements to derive empirical relations to the measured Z, yielding a variety of Z–S relations, depending on the climatological properties of snow at the measurement location.

Improvements upon these situational Z–S relationships can be made if multiparameter radar observations are available. For snow, these have historically proceeded along two pathways following advances in multi-frequency/Doppler and dual-polarization radar technologies. Multi-frequency methods essentially rely upon deviations from Rayleigh scattering to infer a characteristic particle size (e.g., Matrosov et al., 2005; Liao et al., 2016) and, with three frequencies (typically X or Ku, Ka, and W bands), density can also be inferred (Kneifel et al., 2015). For vertically pointing radars, Doppler velocities can also be used to refine the density estimate (Oue et al., 2015; Mason et al., 2018), since fall velocity of a snow particle depends (to first order) on its size and density (Heymsfield and Westbrook, 2010). These methods have been employed primarily towards data collected at a few well-equipped snow observatories such as the Hyytiälä Forestry Research Station in Hyytiälä, Finland (Hari and Kulmala, 2005), the Jülich Observatory for Cloud Evolution (JOYCE) in Germany (Löhnert et al., 2015), and the Department of Energy – Atmospheric Radiation Measurement Program facility in Alaska (de Boer et al., 2018) and applied to airborne and spaceborne radar datasets (Leinonen et al., 2018; Tridon et al., 2019; Chase et al., 2021).

Polarimetric radar measurements have shown value in inferring ongoing ice growth processes, due to the dependence of these measurements on the distributions of particle shapes, orientations, and sizes. In particular, enhancements in the differential reflectivity (Z_DR) and specific differential phase (K_dp) have been linked to the planar crystal growth near −15 ^∘C (e.g., Ryzhkov and Zrnić, 1998; Kennedy and Rutledge, 2011; Andrić et al., 2013; Schrom et al., 2015; Moisseev et al., 2015). Assessing changes in vertical profiles of the polarimetric radar variables also provides information on ice growth processes. Decreases in Z_DR and K_dp towards the ground have been observed with increases in reflectivity towards the ground, indicating growing ice particles becoming more chaotically oriented and more spherical, a result of some combination of aggregation and intense riming (e.g., Bechini et al., 2013; Oue et al., 2015; Ryzhkov et al., 2016; Schrom and Kumjian, 2016; Kumjian and Lombardo, 2017). Some recent efforts have been made to gain quantitative information about the ice particle properties (and thus associated microphysical processes and snowfall rates) from polarimetric radar measurements using empirically determined algorithms (e.g., Bukovčić et al., 2020) and microphysical-model informed parameter estimation (e.g., Schrom et al., 2021). However, there has been limited evaluation of these methods using additional remote sensing (e.g., multi-frequency radar measurements) and in situ observations.

From the extensive literature on multi-frequency and polarimetric radar studies of snow, it is evident that complementary information is contained in these measurements. However, relatively few studies have been performed to assess the information content of multi-frequency, dual-polarization radar measurements of snow, partially because of a lack of radar platforms with these capabilities deployed in locations subject to frequent and high-accumulation snowfall events. The NASA dual-polarization, dual-frequency Doppler radar (D3R) is a premier radar for making such measurements. The D3R was built to provide ground validation measurements for the NASA Global Precipitation Measurement (GPM) mission's dual-frequency precipitation radar (DPR; Chandrasekar et al., 2010; Vega et al., 2014) and operates at Ku and Ka bands using novel solid-state transmitters. The D3R was first deployed in the GPM Cold-season Precipitation Experiment (GCPex; Skofronick-Jackson et al., 2015) and in subsequent GPM ground validation field campaigns. Upgrades to improve sensitivity and range resolution have been implemented (Kumar et al., 2017) since these early campaigns.

Recognizing the potential utility of D3R measurements to provide unique information about microphysics, dynamics, and quantitative precipitation estimation (QPE) in snowstorms, the Korean Meteorological Administration (KMA), organizers of the International Collaborative Experiments for PyeongChang 2018 Olympic and Paralympic winter games (ICE-POP 2018) cooperated with NASA to deploy the D3R radar in Daegwallyeong, South Korea, from November 2017–March 2018. The D3R was part of an extensive network of ground-based remote sensing and in situ instrumentation deployed during ICE-POP 2018 and formed a central observation point for measurements aligned perpendicular to the coastal mountain ranges of eastern South Korea. This measurement strategy was devised to examine the distribution of precipitation from the coast to the mountains in different winter synoptic weather situations and evaluate high-resolution numerical weather prediction in this complex topographic region (Lim et al., 2020).

The objective of this study is to use the data collected by the D3R during ICE-POP 2018 develop a snow retrieval algorithm using realistic scattering models of pristine, aggregate, and rimed snow particles to further our understanding of the complementary nature of the dual-frequency and polarimetric radar measurements and their utility regarding snow microphysical characterization and QPE. The output of this algorithm is intended to aid in identifying microphysical processes during ICE-POP events and provide snow QPE during the deployment. The extensive network of ground instrumentation is leveraged to validate the algorithm output. The manuscript is organized as follows: the observational datasets are listed in Sect. 2, the particle scattering properties we use and the construction of lookup tables from these databases is described in Sect. 3, algorithm mechanics and information content are analyzed in Sect. 4, validation for selected cases is presented in Sect. 5, and the conclusions are given in Sect. 6.

2 Datasets

2.1 D3R

The NASA D3R radar is a polarimetric Doppler weather radar operating at Ku (13.91 GHz) and Ka (35.56 GHz) bands, which utilizes novel design features including aligned antennas, solid-state transceivers, and a digital waveform generator to enable deployment in a wide range of environmental conditions on a mobile trailer platform (Chandrasekar et al., 2010). At both frequencies, the following parameters are measured: reflectivity (Z), differential reflectivity (Z_dr), differential propagation phase (ϕ_dp), co-polar correlation coefficient (ρ_hv), radial velocity (V), and spectrum width (W).

During ICE-POP 2018, the D3R was located on the roof of the Daegwallyeong (DGW) regional weather office (36.677^∘ N, 128.719^∘ E; altitude 789 m m.s.l.). The D3R was configured to measure 150 m range gates out to a maximum range of 39.75 km, combining a short pulse for ranges <3.3 km with a medium pulse for the remaining range gates. This gives a Ku-band sensitivity ranging approximately from −30 to −5 dBZ over the short pulse and from −15 to 5 dBZ over the long pulse (Kumar et al., 2017). The primary scan schedule during snow events conducted a Plan Position Indicator (PPI) scan at 5^∘ elevation followed by Range Height Indicator (RHI) scans at 51, 231, and 330^∘ azimuths, with each set of scans taking 5 min. For this study, we focus on the 231^∘ RHI scans aimed towards the May Hills Supersite (MHS) 2 km downrange, which contained a wealth of ground instrumentation.

We use the D3R data available from the NASA ICE-POP data archive held at the Global Hydrology Resource Center (Petersen et al., 2018). Several snowfall events were observed by D3R during the ICE-POP 2018 campaign, but we selected only three events to analyze in this manuscript, listed in Table 1, for their diversity in synoptic forcing, environmental profiles, and microphysics. Time–height profiles from the Ku-band polarimetric RHIs near the Precipitation Imaging Package (PIP) are shown for these in Fig. 1. The 9 January case was relatively shallow and had increases in reflectivity at the lowest altitudes of the radar beam close to the ground (Fig. 1a). The 28 February and 7–8 March cases were deeper, with reflectivity >20 dBZ extending above 2 km. These cases both had larger enhancements in K_dp, with the 28 February having the largest K_dp and Z_H values of the three cases (Fig. 1b and h). We disregarded events that had mixed-phase precipitation, as the modeling of their scattering properties is less mature than ice particles for our needs, as well as events that did not have both the D3R and data from ground instruments at MHS available.

Table 1 Selected snowfall events during ICE-POP 2018. Synoptic classification is from Kim et al. (2021). SWER indicates snow water equivalent rate.

Figure 1 Time–height plots of vertical profiles of Z_H (first row), Z_DR (second row), K_dp (third row), ρ_hv (fourth row), and DWR (fifth row) for the 9 January 2018 (left column), 28 February 2018 (middle column), and 7–8 March 2018 (right column) cases.

The D3R Z_dr and ϕ_dp were calibrated on an event basis, and the absolution calibration of Ku and Ka Z was established in the early period of the deployment (Chandrasekar et al., 2018). Notwithstanding these calibration efforts, in order for the retrieval algorithm to perform optimally, we examined the data for self-consistency of Z_dr and the dual-wavelength ratio (DWR), defined as Z_Ku−Z_Ka. The expectation is that, at near-zenith angles, Z_dr should be close to zero, and, in the limit of small particles, the DWR should equal the difference in attenuation, which should be small, but positive, depending on the path-integrated attenuation from water vapor, cloud liquid water, and hydrometeors.

We examined time series of the PDFs of these quantities for the events listed in Table 1. In Fig. 2, the top two rows show the time series of the PDF of Ku and Ka Z_dr at elevation angles $>\mathrm{80}{}^{\circ }$ and altitudes above 1.5 km m.s.l. There is clearly some non-stationary behavior in these PDFs, so we applied additional calibration offsets (independently at each frequency) to the Z_dr for each RHI such that the average Z_dr at elevation angles $>\mathrm{80}{}^{\circ }$ is equal to zero. The time series of the calibrated Ku and Ka Z_dr are shown in the third and fourth rows of Fig. 2. In addition to the modification to Z_dr, we also noticed some unusual behavior in the small-particle DWR PDF (fifth row of Fig. 2) during the 28 February case. When Z_Ku<10 dBZ, Skofronick-Jackson et al. (2019) found that Ku–Ka DWR is typically close to zero, but during this case, there is large increase in the DWR, peaking around 06:00 UTC, that is difficult to explain entirely by particle size or attenuation. Accumulation of wet snow on the radome is a possible explanation for this behavior (Venkatachalam Chandrasekhar, personal communication, 2020), but for this study we have chosen to avoid correcting the DWR because it is difficult to have a continuous independent estimate of the Ku–Ka relative calibration from a ground radar.

Figure 2 Time-series histograms of near-zenith Ku Z_dr before (top row) and after (second row) offsets were applied. Similar histograms are presented for Ka-band data in rows three and four. The fifth row displays the time-series histogram of the Ku–Ka DWR where Ku reflectivity is less than 10 dBZ within 10 km of the radar. Colors are indicative of relative occurrence, normalized at each time step.

2.2 PIP

The primary source of validation for the microphysical retrievals in this study is PIP (Pettersen et al., 2020), which takes measurements that can be used to derive quantities including snowfall rate and effective density. The PIP is a video disdrometer made up of a single high-speed camera, continuously recording at 380 frames per second, and a halogen lamp, which is used to backlight the precipitation particles. The camera and lamp are separated by 2 m and the focal plane is located 1.33 m from the camera lens and uses an open sampling volume (i.e., the sampling volume is not enclosed within a box). The images produced are 640×480 pixels with a resolution of 0.1 mm×0.1 mm. The field of view (FOV), including the edge effects, is then 64−D_eq mm $\times \mathrm{48}-D_{\mathrm{eq}}$ mm, where D_eq is the equivalent diameter in millimeters. The depth of field (DOF) also varies with particle size and is expressed as $\mathrm{117}/D_{\mathrm{eq}}$ (in mm). The sampling volume is a multiplication of the FOV, DOF, and the number of frames over a given time period. Considering 100 particles with uniform size of D_eq=1 mm, the sampling volume is 790 m³ for a 1 min observation period. As PIP only uses a single camera, the precipitation particle images are of a projected view of the particle and do not contain any information on the particle dimension along the viewing direction.

PIP determines the characteristics of the precipitation particles using an algorithm written using the National Instruments IMage AQuisition (IMAQ) software package. This algorithm determines the shape of the precipitation particle (i.e., the long and the short dimensions) by fitting an ellipse to the PIP-imaged particle. The IMAQ software package defines the fitted ellipse as the ellipse having both the same area and the same perimeter as the PIP-imaged particle. During our preliminary analysis of the data, we found that the IMAQ-fitted ellipses tended to overestimate the long dimension of the particle and underestimate the short dimension, resulting in an underestimate of the aspect ratio. As such, we have reprocessed the PIP data using an alternative, custom-built ellipse-fitting strategy (Helms et al., 2022). This strategy uses the method implemented in the fit_ellipse program in the Coyote Interactive Data Language (IDL) library (http://www.idlcoyote.com/, last access: 25 February 2022) to perform the actual fit. The fitting is performed on images of particles taken from videos that PIP records for troubleshooting purposes. These videos contain the first 2000 frames that contain precipitation particles within each 10 min period. For periods with fewer than 2000 frames containing precipitation particles, the videos will be shorter than 2000 frames.

The PIP-determined orientation angle is defined as the counterclockwise angle from the positive horizontal axis, where positive is to the right, to the longest dimension of the particle. This results in orientation angles ranging from 0 to 180^∘. In order to combine the ellipse aspect ratio (minor axis length divided by major axis length) and orientation angle information, we have used the natural logarithm of the aspect ratio of the bounding box of the particle. The bounding box is defined as the smallest rectangle that is able to contain the particle and whose edges are either horizontally or vertically oriented.

2.3 Soundings

Radiosonde launches were performed every 3 h during precipitation events from the DGW Regional Weather Office, supplementing the normal 12-hourly observations, and used Meteomodem M10 radiosondes (Meteomodem, 2021). The profiles of temperature and humidity from the nearest-in-time radiosonde were used as input to the D3R algorithm to calculate attenuation from atmospheric gases. These profiles are also used to qualitatively evaluate the retrieval output with respect to locations of well-known thermodynamic importance (e.g., the −15 ^∘C dendritic growth maximum and layers that are supersaturated with respect to ice).

3 Particle scattering properties

Gaining useful information from polarimetric radar measurements requires scattering properties of hydrometeors with preferred mean orientations. We incorporate such scattering properties into the retrieval algorithm described herein, by way of lookup tables (LUTs) that are derived by integrating these scattering properties over a prescribed PSD. Specifically, we use scattering properties for a pristine plate, an aggregate of plates, and graupel (all at a wide range of sizes) from the Atmospheric Radiative Transfer Simulator (ARTS; Eriksson et al., 2018; Brath et al., 2020; Ekelund et al., 2020) database. Scattering properties for each of these particles are available for a discrete set of frequencies (for this study, we use the calculations at 13.4 and 35.6 GHz), temperatures (190, 230, and 270 K), incident angles of the transmitted radiation, scattering angles of the scattered radiation, and, except for graupel (which is assumed to have total random orientation) zenith-relative orientation angles (we refer to this angle hereafter as β). Table 2 lists values of these properties that correspond to the scattering calculations.

Table 2 Properties of the ARTS scattering database calculations used herein.

NA: not available.

The ARTS database aims to provide scattering properties for a set of particles over a large range of frequencies so that applications using both active and passive remote sensing measurements are consistent. The derivation of polarimetric radar variables from this database is described in Appendix A. For each particle type listed in Table 2, these variables are integrated over a PSD using a modified gamma form (e.g., Petty and Huang, 2011):

$$\begin{array}{}\text{(1)}& N\left(D\right)=N_{\mathrm{0}}{D}^{\mathit{\mu }}\exp \left[-{\displaystyle \frac{\mathit{\mu }+\mathrm{4}}{D_{\mathrm{m}}}}D\right],\end{array}$$

where D is the diameter of an equivalent-volume solid ice sphere in mm, D_m is the mass-weighted mean equivalent-volume diameter, and μ is the shape parameter. While this definition of D (and D_m) complicates the comparison to in situ measurements (where maximum dimension is typically used to describe size), it is directly related to particle mass, and, in the Rayleigh limit, radar reflectivity. To reduce the dimensionality of the LUTs while preserving the variables that control the shape of the PSD, only D_m and μ are varied in the construction of the LUT, and N₀(D_m,μ) is a normalized concentration factor such that the ice water content of all PSDs stored in the LUT is 1 g m⁻³:

$$\begin{array}{}\text{(2)}& N_{\mathrm{0}}\left(D_{\mathrm{m}},\mathit{\mu }\right)={\displaystyle \frac{\mathrm{1}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}}{{\int }_{D_{min}}^{D_{max}}\frac{\mathit{\pi }}{\mathrm{6}}{\mathit{\rho }}_{\mathrm{i}}{D}^{\mathrm{3}+\mathit{\mu }}N\exp \left[-\frac{\mathit{\mu }+\mathrm{4}}{D_{\mathrm{m}}}\right]\mathrm{d}D}},\end{array}$$

where ρ_i is the density of solid ice (0.917 g cm⁻³). The distribution is later scaled by a retrieved concentration factor to match the observed reflectivity. For the particle types with preferred orientation, these PSDs are further integrated over the range of β angles to account for the orientation distribution. We choose the von Mises distribution to represent the PDF of beta angles (Table 2). The von Mises distribution is a continuous function that represents the dispersion of variables in circular coordinates, and has been used to represent hydrometeor orientation retrieved from polarimetric radar (Bringi and Chandrasekar, 2001; Melnikov and Straka, 2013). For a zero-mean canting angle, this distribution simplifies to

$$\begin{array}{}\text{(3)}& p\left(\mathit{\beta }\right)={\displaystyle \frac{e^{\mathit{\kappa }\cos \mathit{\beta }}}{G\left(\mathit{\kappa }\right)}},\end{array}$$

where κ is a dispersion parameter and G(κ) is a normalization factor such that the sum of probabilities is equal to 1. When κ=0, the distribution is uniform; as κ increases, the distribution becomes narrower and can be approximated by a normal distribution with standard deviation $\sqrt{\mathrm{1}/\mathit{\kappa }}$. In the construction of the LUTs, we make the simplifying assumption that κ is independent of particle size. While this is almost certainly an oversimplification (observations and Reynolds number analysis suggest that smaller particles of a given aspect ratio should have more randomly distributed canting angles; Klett, 1995), in the retrieval, this is dealt with by allowing the combination of pristine and aggregate PSDs of different κ values, as will be further described in Sect. 4.

Table 3 Dimensions of the lookup tables derived from the particle types selected from the ARTS single-scattering database.

Table 4 Variables stored within the lookup tables.

The LUTs are constructed for each particle species. The LUT dimensions and ranges of the LUT indices are given in Table 3. The variables stored within the LUTs encompass three categories: radar variables, physical variables, and simulated PIP measurements. The radar variables are the single-particle scattering properties necessary to construct the polarimetric radar measurements: the reflectivity factor Z at both horizontal and vertical polarizations, the extinction coefficient k at both polarizations, the specific differential phase K_dp, and the real and imaginary parts of the co-polar conjugate product of scattering amplitudes ($C_{\mathrm{hv}}=S_{\mathrm{hh}}^{*}{S}_{\mathrm{vv}}$). The formulas used to forward model the radar measurements from these quantities are given in Appendix B.

In addition to the radar variables, we simulate several variables (denoted with subscript P) that represent PIP-measured quantities. Because the PIP measures 2-D projections of 3-D particles, assumptions must be made to convert the PIP measurements to physical variables. While several formulas have been proposed to achieve this purpose (e.g., Jiang et al., 2017), we opted instead to use 2-D projections of particles derived from the shapefiles used in the ARTS database to simulate the PIP measurements over the range of D_m, μ, and κ in the LUTs, as this process is less ambiguous than the alternative of attempting to derive the 3-D properties from the 2-D PIP measurements. Moreover, this process ensures an internally consistent comparison of the radar retrieval output with PIP measurements. Most of these quantities depend in some way on the equivalent diameter D_eq measured by the PIP. The simulated PIP measurements include the snowfall water equivalent rate (S), volumetric snowfall rate (S_V), effective density (ρ_P), normalized intercept $N_{\mathrm{P}}^{*}$, area-weighted mean particle volume diameter D_P, mean box aspect ratio ($\stackrel{\mathrm{‾}}{a_{\mathrm{box}}}$), and mean ellipse aspect ratio ($\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$). The final two quantities are weighted by the projected area. The snowfall water equivalent rate is calculated as

$$\begin{array}{}\text{(4)}& S=\mathrm{3.6}\underset{\mathit{\beta }}{\int }p\left(\mathit{\beta }\right)\underset{D_{min}}{\overset{D_{max}}{\int }}N\left(D\right)m\left(D\right)V_{\mathrm{t}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta },\end{array}$$

where m(D) is the mass of the particle in grams, V_t(D,β) is the terminal velocity in m s⁻¹, and S is in units of mm h⁻¹. The terminal velocity is calculated following the process given in Heymsfield and Westbrook (2010), where the area ratio is simulated from the ARTS shapefiles (and thus depends on both D and β), and using values for the density and viscosity of air typical at the MHS location during ICE-POP 2018 snow events (−5 ^∘C, 90 %RH, 925 hPa). The volumetric snowfall rate is calculated using the same terminal velocity but, instead of mass, using the volume of the particle derived from the PIP-measured diameter:

$$\begin{array}{}\text{(5)}& S_{\mathrm{V}}=\mathrm{3.6}\underset{\mathit{\beta }}{\int }p\left(\mathit{\beta }\right)\underset{D_{min}}{\overset{D_{max}}{\int }}N\left(D\right){\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{D}_{\mathrm{eq}}^{\mathrm{3}}(D,\mathit{\beta })V_{\mathrm{t}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }.\end{array}$$

The PIP-derived density (ρ_P) is defined as the ratio of the liquid-equivalent snow rate to the volumetric snowfall rate multiplied by the density of liquid water ρ_w (Tiira et al., 2016):

$$\begin{array}{}\text{(6)}& {\mathit{\rho }}_{\mathrm{P}}={\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{S}{S_{\mathrm{P}}}}.\end{array}$$

The normalized intercept $N_{\mathrm{P}}^{*}$ is adapted from the definition given for $N_{\mathrm{0},\mathrm{23}}^{*}$ in Field et al. (2005):

$$\begin{array}{}\text{(7)}& N_{\mathrm{P}}^{*}={\displaystyle \frac{{\left[{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)d_{\max \mathrm{2}\mathrm{D}}^{\mathrm{2}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }\right]}^{\mathrm{4}}}{{\left[{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)d_{\max \mathrm{2}\mathrm{D}}^{\mathrm{3}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }\right]}^{\mathrm{3}}}},\end{array}$$

where d_max2D(D,β) is the average 2-D projected maximum particle dimension. This parameter is useful for providing temperature-dependent constraints on the PSD as will be shown in Sect. 4.

D_P is the ratio of the fourth to third moments of the PSD in terms of D_eq:

$$\begin{array}{}\text{(8)}& D_{\mathrm{P}}={\displaystyle \frac{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{4}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{3}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}}.\end{array}$$

In order to evaluate the capability of the D3R polarimetric measurements to retrieve a bulk measurement of the aspect ratio, we defined the area-weighted mean ellipse aspect ratio ($\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$):

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}={\displaystyle \frac{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{2}}(D,\mathit{\beta })\frac{b}{a}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{2}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}},\end{array}$$

where $\frac{b}{a}(D,\mathit{\beta })$ is the ratio of the short to long axis of the ellipse fitted to the simulated PIP image. In order to compare our retrievals to a variable that is influenced by the canting angle distribution parameter κ, we defined the area-weighted mean box aspect ratio ($\stackrel{\mathrm{‾}}{a_{\mathrm{box}}}$):

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{a_{\mathrm{box}}}={\displaystyle \frac{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{2}}(D,\mathit{\beta })\frac{y}{x}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}{{\int }_{\mathit{\beta }}p\left(\mathit{\beta }\right){\int }_{D_{min}}^{D_{max}}N\left(D\right)D_{\mathrm{eq}}^{\mathrm{2}}(D,\mathit{\beta })\mathrm{d}D\mathrm{d}\mathit{\beta }}},\end{array}$$

where $\frac{y}{x}(D,\mathit{\beta })$ is the ratio of the vertical to horizontal dimensions of the bounding box of the simulated PIP image.

While it is not feasible to illustrate every dimension of the LUTs, a few examples are provided to preview the expected sensitivity of the D3R measurements to microphysical parameters. In Fig. 3, the Ku–Ka dual-wavelength ratio is plotted as a function of both the volume-equivalent D_m and PIP-measured D_P for each of the three species. The DWR is nearly identical for the pristine plates and plate aggregates at D_m<1 mm and D_P<3 mm but reaches an upper limit of about 5 dBZ for the single plates while increasing to nearly 10 dBZ for the aggregates. Large Ku–Ka DWRs, sometimes exceeding 10 dBZ, have been observed coincident with large aggregates (20 mm projected diameter) in dendritic growth regimes during OLYMPEX (Chase et al., 2018), so it is important that the LUTs capture this DWR range. The DWR of graupel is lower than the unrimed plates and aggregates for a given volume-equivalent size but is similar or greater when viewed with respect to PIP-measured size based on the projected area. This suggests that DWR alone is not sufficient to determine the PSD mean particle size and density. Such information can be provided by the polarimetric measurements. For this study, we assume graupel to be randomly oriented and thus have zero contribution to Z_dr and K_dp, which is a reasonable assumption for dry graupel (Kumjian, 2013). Meanwhile, the contributions to Z_dr and K_dp from plates and aggregates primarily depend on D_m and κ (Fig. 4). As expected, both Z_dr and K_dp increase with κ as the orientation distribution becomes less dispersive. The Z_dr rapidly decreases with aggregate size, as does K_dp (though less rapidly). Meanwhile, the Z_dr of the individual plates does not depend on D_m, since these plates have a fixed aspect ratio at sizes larger than about 0.2 mm (Brath et al., 2020). The K_dp of the individual plates also does not depend much on size, although some resonance effects lead to a small decrease in K_dp at larger size. From Figs. 3 and 4, it is clear that there is complementary information in the dual-frequency and polarimetric variables to discern particle size and species, which will be demonstrated in the next section.

Figure 3 Ku–Ka dual-wavelength ratio for the three ice species listed in Table 2 as a function of the melted (volume equivalent) D_m (a) and D_P (b). The families of curves represent different values for the shape parameter μ ranging from −1 to 5, with lighter shades indicating higher values and less dispersion of the PSD. The different ranges of D_P are a consequence of the different size–density relationships for each species.

Figure 4 Polarimetric variables Z_dr (a, b) and K_dp (c, d) as a function of mean particle size D_m and orientation distribution parameter κ. Panels (a, c) represent aggregates, and (b, d) individual plates. All values are for horizontal incidence at the Ku band at 270 K.

4 Algorithm description

Optimal estimation (OE; Rodgers, 2000) is a form of Bayesian inversion that assumes Gaussian error statistics and accommodates moderately nonlinear forward models. OE has been used for single-frequency (L'Ecuyer and Stephens, 2002; Munchak and Kummerow, 2011) and multi-frequency (Grecu et al., 2011; Mason et al., 2018) radar precipitation retrievals. It is applied here to the multi-frequency, polarimetric observations provided by the D3R. In this section, the OE components (state vector, observation vector, and covariance matrices) are defined and the approach is illustrated with an example retrieval along a single radar ray.

4.1 Optimal estimation setup

The cost function that is minimized by optimal estimation is

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{\Phi }=& \left({\mathit{Y}}_{\mathit{s}\mathit{i}\mathit{m}}\right(\mathit{X})-{\mathit{Y}}_{\mathit{o}\mathit{b}\mathit{s}}{)}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{y}}^{-\mathbf{1}}({\mathit{Y}}_{\mathit{s}\mathit{i}\mathit{m}}\left(\mathit{X}\right)-{\mathit{Y}}_{\mathit{o}\mathit{b}\mathit{s}})\\ & +{\left(\mathit{X}-{\mathit{X}}_{\mathit{a}}\right)}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{a}}^{-\mathbf{1}}\left(\mathit{X}-{\mathit{X}}_{\mathit{a}}\right),\end{array}\end{array}$$

where Y_obs is the observation vector and Y_sim(X) its forward-modeled counterpart, S_y is the measurement and forward model error covariance matrix, X is the state vector and X_a its prior, and S_a is the state covariance matrix. For the D3R retrieval of snow microphysical properties, we define Y_obs to contain one or more of the following series of dual-frequency or dual-polarization observations at each range gate along a radial:

$$\begin{array}{}\text{(12)}& {\mathit{Y}}_{\mathit{o}\mathit{b}\mathit{s}}=\left[\begin{array}{ccc}{\text{DWR}}^{\mathrm{1}}& \mathrm{\cdots }& {\text{DWR}}^{n\mathrm{DWR}}\\ Z_{\mathrm{drKu}}^{\mathrm{1}}& \mathrm{\cdots }& Z_{\mathrm{drKu}}^{n\mathrm{ZdrKu}}\\ Z_{\mathrm{drKa}}^{\mathrm{1}}& \mathrm{\cdots }& Z_{\mathrm{drKa}}^{n\mathrm{ZdrKa}}\\ {\mathit{\varphi }}_{\mathrm{dpKu}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\varphi }}_{\mathrm{dpKu}}^{n\mathit{\varphi }\mathrm{dpKu}}\\ {\mathit{\varphi }}_{\mathrm{dpKa}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\varphi }}_{\mathrm{dpKa}}^{n\mathit{\varphi }\mathrm{dpKa}}\end{array}\right],\end{array}$$

where each type of observation is filtered for ground clutter and only considered if the signal-to-noise ratio exceeds 1. Thus, the number of observations of each type may differ. Note that ϕ_dp instead of K_dp was chosen because ϕ_dp is the more direct measurement, and additional assumptions (with potentially non-Gaussian errors) are required to derive K_dp from noisy ϕ_dp measurements.

The state vector consists of quantities that describe the PSD, relative contribution of each species, as well as other quantities known to affect the Ku- and Ka-band polarimetric radar measurements:

$$\begin{array}{}\text{(13)}& \mathit{X}=\left[\begin{array}{ccc}\mathit{\delta }{N}_{\mathrm{P}}^{*\mathrm{1}}& \mathrm{\cdots }& \mathit{\delta }{N}_{\mathrm{P}}^{*n\mathrm{NP}}\\ f_{\mathrm{p}}^{\mathrm{1}}& \mathrm{\cdots }& f_{\mathrm{p}}^{n\mathrm{fp}}\\ f_{\mathrm{r}}^{\mathrm{1}}& \mathrm{\cdots }& f_{\mathrm{r}}^{n\mathrm{fr}}\\ D_{\mathrm{p}}^{\mathrm{1}}& \mathrm{\cdots }& D_{\mathrm{p}}^{n\mathrm{Dp}}\\ {\mathit{\kappa }}_{\mathrm{p}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\kappa }}_{\mathrm{p}}^{n\mathit{\kappa }\mathrm{p}}\\ {\mathit{\kappa }}_{\mathrm{a}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\kappa }}_{\mathrm{a}}^{n\mathit{\kappa }\mathrm{a}}\\ {\mathit{\mu }}_{\mathrm{p}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\mu }}_{\mathrm{p}}^{n\mathit{\mu }\mathrm{p}}\\ {\mathit{\mu }}_{\mathrm{a}}^{\mathrm{1}}& \mathrm{\cdots }& {\mathit{\mu }}_{\mathrm{a}}^{n\mathit{\mu }\mathrm{a}}\\ c^{\mathrm{1}}& \mathrm{\cdots }& c^{n\mathrm{c}}\end{array}\right],\end{array}$$

where each quantity is defined at nodes (indicated by the superscript in Eq. 13) which may be arbitrarily placed (for this study, nodes are spaced 600 m horizontally and 250 m vertically). Each of these quantities, their priors, and ranges are defined in Table 5. Following Grecu et al. (2011), the measured Ku-band Z_hh (which is not in Y_obs) is used as direct input to the forward model, and from these measurements and the quantities in the state vector X, the measurements in Y_obs are simulated. A detailed description of this forward model is provided in Appendix B.

Field et al. (2005)

Table 5 Names, definitions, a priori value, standard deviation, and allowed ranges of quantities in the state vector (Eq. 13). Note that some quantities are retrieved in logarithmic space to accommodate lognormal distributions and/or to linearize the forward model. SD: standard deviation.

The observation error covariance matrix S_y must accurately describe the error characteristics of the measurements and forward model. Values that are unrealistically low can lead to overfitting, whereas overly conservative (large) error estimates can lead to underutilization of the information contained within the measurements. In this study, we consider the diagonal elements of S_y to be the sum of the measurement error and forward model error components, which are given in Table 6 for each measurement in S_y. The measurement errors are obtained from the gate-to-gate variance over many homogeneous scenes at low (<10^∘) elevation angles, where the assumption is that the true change in the measured quantity is small compared to the measurement noise. The measurement error is assumed to be uncorrelated in space and between variables. The forward model error is quantified by assessing the variance in the measurements for alternate aggregate particles of the same size as the ARTS aggregates. These alternate models include dendrites and columns, with different assumptions about the aggregation process (Schrom et al., 2021). Although these forward model errors are correlated between variables, analysis of the covariance matrices showed the correlation between DWR and Z_dr error to be insignificant. While Z_dr−K_dp error correlation is larger, the ϕ_dp error is dominated by upstream propagation error which is assumed to be uncorrelated to the Z_dr forward model error at a given range gate. Therefore, for this study, the off-diagonal elements of S_y are set to zero, although this is a choice that could be refined in future implementations of the retrieval.

Table 6 Measurement and forward model components of the error covariance matrix S_y for each measurement, expressed as standard deviations.

4.2 Example ray

The optimal estimation process along a ray of radar data is illustrated in Fig. 5, which shows the observed (Y_obs) and simulated (Y_sim) measurements, and Fig. 6, which shows the various retrieval parameters in X as well as derived quantities of ice water content and D_m for each iteration until convergence. This ray is characterized by DWR peaks at 0–5 and 23–28 km, which reach 6 dB, dropping to 1–3 dB elsewhere. The Z_dr has several peaks, the most significant of which reach over 2 dB at 20 and 30 km range, bracketing the DWR peak. The ϕ_dp increases most rapidly in the 10–20 km range with smaller rates of increase elsewhere, reaching 15^∘ at the Ku band and 50^∘ at the Ka band. The initial profiles of the retrieval parameters $N_{\mathrm{P}}^{*}$, f_r, f_p, κ_p, κ_a, μ_p, μ_a, and the derived ice and cloud liquid water content and D_m of the aggregate and pristine PSDs are shown in the lightest shaded lines in Fig. 6. The iterative adjustments guided by the Jacobian matrix respond to the initial differences between Y_obs and Y_sim:

• $N_{\mathrm{P}}^{*}$ decreases slightly near the DWR peaks and increases elsewhere.

• f_p decreases below 50 % in the DWR peak region but increases above 60 % in the 10–20 km range and again near 30 km, corresponding to the Z_dr peaks and steepest increase in ϕ_dp.

• f_r is generally below 40 %, with minima near the Z_dr peaks.

• κ_p is generally lower than κ_a but exhibits peaks corresponding to the Z_dr peaks, whereas κ_a does not vary as much with range.

• μ_a and μ_p increase slightly from their initial values in most regions, although μ_a dips below zero near the DWR peaks (lower μ corresponds to a longer tail of the PSD at large sizes and can result in higher DWRs, all else being equal).

• No significant cloud liquid water is detected or needed to explain the DWR observations. In fact, the near-zero DWR at the far ranges implies that there is little differential attenuation, and the observed near-field DWR can be attributed to particle size effects.

• D_m of the pristine population (controlled by the retrieval parameter D_p) generally increases as a proportion of the D_m of the aggregates in order to balance pristine particle contributions to Z_dr and ϕ_dp.

Figure 5 D3R observations of reflectivity (a) and DWR (d), Z_dr (b, e), and ϕ_dp (c, f) along the 231^∘ azimuth, 6^∘ elevation radial at 02:17 UTC on 8 March 2018. For the DWR, Z_dr, and ϕ_dp, the D3R observations are indicated by the solid black lines and the simulated measurements for each iteration are shown in progressively more opaque blue lines, with the final iteration plotted in red.

The most significant impact of these parameter changes is to increase D_m and reduce ice water content in the DWR peak regions, with the opposite changes elsewhere. The final retrieved state is in good agreement with the Z_dr and ϕ_dp observations at both frequencies. There is a residual high bias in the DWR, especially at the range beyond 30 km. Examining the DWR histogram for this case in Fig. 2 reveals a mean near or below zero, which may indicate a low bias in the relative calibration of the Ku to Ka reflectivity, and possible high bias in the ice water content retrieval.

Figure 6 Retrieved $N_{\mathrm{P}}^{*}$ (a), f_r and f_p (b), κ_p and κ_a (c), μ_p and μ_a (d), water content (e), and D_m (f) from the D3R observations of Fig. 5. As in Fig. 5, the retrieved measurements for each iteration are shown in progressively more opaque lines.

As a consistency check, we show that the retrieval algorithm effectively reproduces the spatial distribution of the radar variables for an RHI (Fig. 7). The retrieval that uses only the DWR as input is able to reproduce the reflectivity and dual-wavelength ratio measurements but overestimates the Φ_dp at both wavelengths and poorly represents the spatial structure of the Z_DR. In contrast, the retrieval using only the Ku-band polarimetric measurements produces simulated Φ_dp and Z_DR that correspond well to the measurements but dual-wavelength ratio simulations that fail to capture the observed dual-wavelength ratio structure in the RHI. The simulated radar variables associated with the retrieval that incorporates all of the radar observations shows the best consistency with the observations, suggesting that the polarimetric and DWR measurements contain independent and complementary information.

Figure 7 Observed and simulated Z_H (first column), Ku-band Z_DR (second column), DWR (third column), Ku-band Φ_dp (fourth column), and Ka-band Φ_dp (fifth column). The observed radar variables are shown on the first row, the simulated radar variables from the Ku-band reflectivity retrieval are shown in the second row, the simulated radar variables from the Ku-band polarimetric retrieval are shown in the third row, the simulated radar variables from the Ku-band and Ka-band reflectivity retrieval are shown in the fourth row, and the simulated radar variables from the full-observation retrieval are shown in the fifth row. The observations are from an RHI taken at 13:00 UTC on 28 February at an azimuth angle of 231^∘.

4.3 Information content analysis

Some further insight into the adjustments made to the parameters can be gained by examining the Jacobian matrix K of partial derivatives of each element of Y_sim with respect to each element of X.

While the Jacobian is state dependent, the sign and relative magnitude of each element from the example ray shown in the previous section are illustrative of the general sensitivity of the forward model to the retrieval parameters. The Jacobian matrix is composed of blocks that are either diagonal or triangular, depending on whether the parameter has a significant downrange propagation effect on an observation. The effect of modifying the parameters in X can be explained by considering the change to the PSD at a fixed Ku-band reflectivity:

• Increasing $N_{\mathrm{P}}^{*}$ results in a smaller mean particle size and higher ice water content. This decreases the DWR and increases ϕ_dp downrange, while the effect on Z_dr is relatively small and state dependent.

• Increasing f_p increases Z_dr and ϕ_dp downrange as a result of the increasing contribution of the pristine particles to the PSD, with little effect on DWR.

• Increasing D_p also increases Z_dr (and decreases DWR) as the pristine particles become larger and contributes more to reflectivity but decreases ϕ_dp downrange as the ice water content is reduced.

• Increasing κ_a and κ_p increases the Z_dr and downrange ϕ_dp, with the κ_p having a much larger magnitude effect.

• Increasing μ_a reduces the DWR as the large-particle tail of the PSD is truncated. There is almost no impact of changing μ_p on the simulated observations.

• Increasing cloud liquid results in an increase in the downrange DWR due to increased differential attenuation but has no impact on the polarimetric measurements.

• Increasing f_r reduces the DWR and all of the polarimetric measurements, as the rimed particles are assumed to have random orientation.

It is noteworthy that the sign of the observation response to perturbations varies in different ways for the different parameters in X. This is an indication that this optimal estimation retrieval, as we have formulated it, is well determined and is also a requirement for reducing ambiguity, or cross talk, in the retrieved state. The information content and cross talk among parameters can also be evaluated by examining the averaging kernel matrix of the retrieved state. The averaging kernel provides a measure of influence of the observations on the retrieved state and is defined as

$$\begin{array}{}\text{(14)}& \mathbf{A}=({{\mathbf{S}}_{\mathbf{a}}}^{-\mathrm{1}}+{\mathbf{K}}^{\mathbf{T}}{{\mathbf{S}}_{\mathbf{y}}}^{-\mathrm{1}}\mathbf{K}{)}^{-\mathrm{1}}({\mathbf{K}}^{\mathbf{T}}{{\mathbf{S}}_{\mathbf{y}}}^{-\mathrm{1}}\mathbf{K}).\end{array}$$

Figure 8 Jacobian matrix of partial derivatives of the simulated D3R measurements to the retrieval parameters. The partial derivatives have been normalized element wise by the square root corresponding diagonal element of S_y – i.e., the expected standard deviation of that measurement's error.

Values close to 1 indicate strong influence of measurements, and values close to 0 indicate that the retrieval is heavily influenced by the prior. In Fig. 9, the median, 10 %, and 90 % quantiles of the averaging kernel are plotted. These statistics were derived from many retrieved states spanning the cases we examined in Table 1. The parameter with the highest information content is the normalized intercept ($N_{\mathrm{P}}^{*}$), which both the DWR and ϕ_dp are highly sensitive to (Fig. 8). The two parameters that describe the pristine component of the PSD (f_p and D_p) also have a similar median and 90 % quantile values to $N_{\mathrm{P}}^{*}$, due to their sensitivity to the polarimetric parameters, but a lower 10 % quantile, likely originating from situation where the f_r is high and there is little sensitivity of the polarimetric variables to f_p and D_p. Similarly, κ_p has a large range between the 10 % and 90 % quantiles, indicating that the information content of this parameter is state dependent, and high values of f_p and D_p are required to maximize the sensitivity of the polarimetric observations to this variable. Another parameter with a wide range of averaging kernel diagonal values is f_r, which requires low f_p and D_p values for it to be the primary driver of the polarimetric variables. Some of this state dependence of the information content is also reflected in the off-diagonal values, which indicate significant cross talk, i.e., a strong correlation in the posterior state vector, between the following groups of variables: $N_{\mathrm{P}}^{*}$, μ_a, f_r; and f_p, D_p, κ_p. These groups have similar Jacobians making it difficult to determine them independently. Finally, it is notable that a few of the variables (κ_a, μ_p, and cloud liquid) have very low averaging kernel values, indicating that the observations are not particularly sensitive to them. This is not surprising, since the aggregates do not show much Z_dr or K_dp sensitivity to κ except at the smallest sizes (Fig. 4), and the shape parameter (μ) of the pristine PSD is not going to influence the reflectivity observations (DWR and Z_dr) because the shape of the large-particle tail is only important to these observations if D_p is close to its upper limit of 1. The sensitivity of Z_dr and K_dp to κ does depend on the aggregate shapes. The ARTS large plate aggregate used herein may not represent cases of more horizontally exaggerated aggregates where the assumed orientation will have a larger impact on the simulated polarimetric variables. The low averaging kernel values for cloud liquid are reflective of the result that it was rarely retrieved in significant quantities, and since it is treated logarithmically, only large values will induce large changes to the DWR. However, there were several RHI scans, particularly on 28 February, where some cloud liquid was necessary to explain high DWR values that could not be achieved by differential scattering alone.

Figure 9 Quantiles of the averaging kernel matrix diagonal elements and off-diagonal elements representing different parameters at the same range gate.

5 Validation

The primary tool used for validation of the D3R retrievals is the PIP located at the MHS location. The D3R retrievals were matched to the PIP by averaging the retrieved quantities in a 600 m wide by 500 m tall box centered above MHS. The lower altitude limit of this box was placed 250 m above the surface to avoid ground clutter contamination along the radials used for averaging. To evaluate the impact of the dual-frequency and dual-polarization measurements, four retrieval experiments were conducted:

• Ku-only: a single-frequency retrieval, equivalent to prescribing a temperature-dependent Z−S relationship;

• DWR-only: a dual-frequency retrieval without polarimetric information; similar information content to the GPM dual-frequency precipitation radar (DPR);

• Ku-pol: a single-frequency polarimetric retrieval using Z_dr and ϕ_dp at the Ku band only;

• All-obs: a dual-frequency polarimetric retrieval using DWR and the Z_dr and ϕ_dp at both frequencies.

To assess the quality of the PIP–D3R matchups, the Ku-band reflectivity was calculated directly from the PIP-derived PSDs using Mie theory (spherical particles) with the PIP-derived particle densities (Tokay et al., 2022). Although there will be some departures from Mie theory for non-spherical particles, at the Ku band these are relatively small (Kuo et al., 2016). The larger contribution to error is the various assumptions required to derive particle density from the PIP size and fall speed (Tokay et al., 2022). Using an ensemble approach to these assumptions, a range of reflectivities was obtained and compared to the D3R-observed reflectivity in the averaging box (top row of Fig. 10).

Some different tendencies are noted for each event. The 9 January case had a consistently higher PIP-derived reflectivity than the D3R measurement, especially during the middle hours of the event. This was a low-echo-top cold low case, and D3R-derived time–height profiles (Fig. 1) indicate significant echo growth (perhaps due to aggregation) at low levels, which may have continued in the clutter region. The 28 February case exhibits a similar bias between the PIP-derived and D3R-measured reflectivity at the Ku band in the first 6 h, after which the D3R reflectivity comes within the lower range of PIP estimates. This was a deeper, warm low case, and the D3R profiles do not indicate any substantial reflectivity increase towards the surface. Instead, it is suspected that excess attenuation from wet snow (the wet bulb temperature was above 0 ^∘C until 04:00 UTC and the air temperature was above 0 ^∘C until 08:00 UTC) accumulated on the radome is responsible for these differences (note also the sharp increase in DWR during the same time period attributed to this factor in Sect. 2). The 7–8 March case exhibited the best agreement between PIP-derived and D3R-measured reflectivity throughout the event. This was also a warm low case with deep echo tops but colder wet bulb temperatures (between −4 and −2 ^∘C during the event).

The validation statistics presented in this section are derived from matchups of the D3R-derived and PIP-measured quantities listed in Table 4. Because the PIP measurements are taken every minute, whereas the D3R RHI scans were conducted every 5–6 min, an optimal lag was found by maximizing the correlation between the lagged D3R-measured and PIP-derived Ku-band reflectivity time series. This lag time was between 2–7 min, depending on the event, which is consistent with a fall speed slightly greater than 1 m s⁻¹ to cover the distance from the center of the averaging box to the surface. Only retrievals where the D3R-measured Ku band reflectivity was within the range of PIP estimates were considered in the calculation of statistics to ensure that the PIP measurements are reasonably representative of the D3R observations.

Figure 10 Time series of D3R-observed and PIP-simulated Ku-band reflectivity (top row), D3R-observed DWR (second row), D3R-observed Z_dr (third row), and D3R-derived K_dp for each event. The PIP simulations used Mie calculations for spherical particles with the PIP-derived effective density. Because these calculations involve a variety of assumptions to derive density from particle shape and fall speed (Tokay et al., 2022), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.

5.1 Snowfall rate and water equivalent

The time series of snow water equivalent rate (S) for each event are shown in Fig. 11. The D3R substantially underestimates snowfall throughout the 9 January and 28 February events, which is not surprising since the PIP-derived reflectivity significantly exceeds the D3R measurement for reasons discussed previously in this section. It is worth noting, however, that the DWR-only and All-obs retrievals are in good agreement with the Pluvio-measured S accumulation on 9 January, which would be consistent with aggregation processes in the lowest levels that increase reflectivity but do not increase S. The 7–8 March event shows better agreement with the PIP measurements. In this case, which had the best reflectivity match, the DWR-only retrieval overestimates S, whereas the Ku-only and Ku-pol retrievals underestimate S relative to the PIP. The retrieval using all of the observations is the closest match and the accumulated S falls within the range of PIP estimates for most of the event. We note from the DWR time series for this event in Fig. 10 that the DWR is just above 0 dB for much of this event which implies a smaller mean particle size and higher S for a given reflectivity. Meanwhile, the Ku-pol retrieval gives a slight reduction in S from the Ku-only retrieval, which is already biased low for this event. Compared to the Pluvio, the All-obs retrieval is biased high; however, after correcting the Pluvio data for wind (Milewska et al., 2019), this retrieval is in better agreement. However, these same wind corrections bias the 9 January event higher than the retrievals and in better agreement with the PIP.

Error statistics for volumetric snowfall rate (S_V) and S from all events, filtered for times when the D3R-measured reflectivity was within the PIP-estimated range, are presented in Table 7. The bias is the overall fraction difference between the accumulated PIP- and D3R-derived amounts. All methods underestimate the amounts, with the DWR-only method providing the closest match for both S_V and S. However, this appears to be the result of compensating biases on the 28 February and 7–8 March events, and the mean absolute error (MAE) is highest for this method. The Ku-pol method provides the best correlation for S_V and S even though it has the largest magnitude bias of all the experiments. The high correlation coefficients, particularly for the multi-parameter methods, are indicative of a response of the radar measurements to the microphysical properties that determine snowfall rate, but the large biases indicate either a calibration bias in the observations, biased forward model (i.e., unrepresentative scattering properties), or both.

Figure 11 Time series of D3R- and PIP-derived snow water equivalent rate (a, b, c) and accumulation (d, e, f) for each event. Because the PIP calculations involve a variety of assumptions to derive particle mass from shape and fall speed (Tokay et al., 2022), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.

Table 7 Mean bias, MAE, and correlation coefficient (r) of PIP-derived versus D3R-estimated volumetric snowfall rate (S_V) and snow water equivalent rate (S) for each of the four retrieval experiments, compiled over all three events when the D3R-observed Ku-band reflectivity was within the PIP-simulated Ku-band reflectivity range. Because of the wide dynamic range of snowfall rates, these quantities are expressed as percentages relative to the mean PIP estimate.

5.2 Mean particle size and density

The D3R retrieval provides an estimate of the PIP-measured area-weighted mean particle volume diameter (D_P) that is consistent with the particle shapes used to generate the LUTs and can be compared directly to the PIP measurement. The time series of this parameter is shown for each event in the top row of Fig. 12. On 9 January 2018, all of the radar retrievals were biased low with respect to the PIP, which is consistent with the reflectivity bias we noted on this day. On 28 February, the PIP measured a rapid increase of D_P to a maximum at 06:00 UTC, followed by a decrease to a minimum around 09:00 UTC and another maximum at 15:00 UTC. None of the retrievals did a particularly good job of capturing the peak at 06:00 UTC, but all methods except the Ku-only retrieval captured the increase toward the second maximum, although again, the peak values were underestimated. On 7–8 March, both methods that used the DWR (DWR-only and All-obs) captured the temporal variability of D_P quite well but were biased low; this is consistent with a suspected low bias in the DWR for this event (mean cloud-top DWRs were slightly below zero; see Fig. 2).

Figure 12 Time series of D3R- and PIP-derived mean volume-weighted particle diameter (a, b, c) and effective density (d, e, f) for each event. Because the PIP calculations involve a variety of assumptions to derive particle density from shape and fall speed (Tokay et al., 2022), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.

The statistics of the D_P retrievals are presented in Table 8 and, as with the snowfall rate statistics, only consider observations where the D3R-observed reflectivity was within the range of calculated PIP values. All of the retrievals are biased low with respect to the PIP, with the Ku-pol retrieval coming the closest. This is counterintuitive, since this retrieval does not consider the DWR, which is the measurement most sensitive to D_P. However, the suspected low bias of the DWR on 7–8 March, which dominates the statistics due to its close match to the observed Ku reflectivity, contributes heavily here. The MAE is only slightly lower for the Ku-pol method than the All-obs method, despite the much more significant bias. Meanwhile, the correlation coefficient is largest for the two methods that incorporate the DWR, suggesting that the DWR is indeed informative regarding the particle size; this underscores the need for DWR to be carefully calibrated to avoid significant bias. Combined use of polarimetric and DWR information seems to at least partially alleviate these biases.

The retrieved effective particle density ρ_P, defined as the snow water equivalent rate divided by the volumetric snowfall rate, can also be directly compared to the PIP measurement (note that this is different from the effective density defined by Pettersen et al., 2020, which is the ratio of fall speeds of an observed snow particle and raindrop of the same equivalent diameter). To first order, ρ_P is inversely proportional to D_P because the unrimed aggregates' density decreases with size. However, changes in f_r and f_p also affect ρ_P, so it is worthwhile to evaluate these retrievals as well. The bottom row of Fig. 12 shows the time series of ρ_P for each event. In general, we find that when the retrieved D_P is biased high with respect to the PIP observation, ρ_P is biased low, and vice versa. It is interesting that on 28 February, all of the D3R methods are in tight agreement, whereas on 7–8 March, the DWR-only and All-obs retrievals are higher than the Ku-only and Ku-pol retrievals, due to the shift toward smaller (denser particles) inferred from the low DWR on this event. In this case, the Ku-only and Ku-pol methods are less biased, but the variability is better represented by the DWR-using retrievals. This is again borne out in the statistics (Table 8), where the Ku-only and Ku-pol methods have the lowest magnitude bias and MAE, and unlike the D_P, the Ku-only method has the best correlation (although the All-obs retrieval is the next best). This suggests that the size–density relationship in the scattering models we chose may not be representative of the particles observed during ICE-POP, or that more information (e.g., W-band reflectivity) is needed to constrain the density (e.g., Kneifel et al., 2015).

Table 8 Mean bias, absolute error (MAE), and correlation coefficient (r) of PIP-derived versus D3R-estimated mean volume-weighted particle diameter (D_vP) and PIP effective density (${\mathit{\rho }}_{\mathrm{P}}=S/S_{\mathrm{V}}$) for each of the four retrieval experiments.

5.3 Bulk particle orientation and aspect ratio

The use of polarimetric measurements in the Ku-pol and All-obs retrievals has the most impact on the retrieved pristine fraction, ratio of pristine to aggregate mean particle size, and pristine population orientation distribution. All of these parameters combine to influence the area-weighted ellipse ($\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$) and box ($\stackrel{\mathrm{‾}}{a_{\mathrm{box}}}$) aspect ratios. In Fig. 13, the influence of the Z_dr and K_dp measurements (Fig. 10) can be observed in these two retrievals, whereas the Ku-only and DWR-only experiments did not appreciably change these parameters. The high Z_dr values before 16:00 UTC on 9 January result in low aspect ratios at that time. There is not much of a trend in Z_dr on 28 February, due to the constant presence of large aggregates which dominate the Z_dr measurement, but there is a notable peak in K_dp around 09:00 UTC which corresponds to the minimum aspect ratio for this event. The 7–8 March event had the least variable polarimetric signatures, but a short peak in Z_dr around 18:00 UTC on 7 March corresponds to a drop in aspect ratio at that time.

Figure 13 Time series of D3R- and PIP-derived area-weighted ellipse (a, b, c) and box (d, e, f) aspect ratio for each event.

The comparison of the retrieved aspect ratios to those derived from the PIP is inconsistent from event to event. There appears to be little variability in the PIP time series on 9 January for either measure of aspect ratio. There is a steady increase in $\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$ between 09:00 and 12:00 UTC on 28 February, matching the retrieved behavior, although the PIP range is considerably smaller than the retrieved range. The 7–8 March event does not show any significant trend in the PIP-measured aspect ratios, consistent with the low Z_dr variability during this event. The aspect ratio error statistics are given in Table 9. The Ku-pol and All-obs methods provide the lowest MAE and least bias; however, the best correlation comes from the Ku-only method (for $\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$) and the DWR-only method (for $\stackrel{\mathrm{‾}}{a_{\mathrm{box}}}$). This is a surprising result, especially in light of the very small range of these retrieved values in Fig. 13. In these methods, the very limited aspect ratio variability is entirely driven by changes in PSD, with smaller mean particle sizes being associated with smaller aspect ratios. The polarimetric measurements add significant variation to the retrieved aspect ratio, but the correlation statistics suggest that the relationship between these measurement and PIP-derived aspect ratio and orientation is tenuous at best. Helms et al. (2022) provide further information on the algorithms used to derive aspect ratio and orientation from the PIP; noting that motion blurring and compression lead to artificially high aspect ratios, particularly for small particles, and quantization artifacts that lead to a maximum aspect ratio of 0.6 for particles less than 0.5 mm in diameter. In future deployments, we expect to measure these properties more precisely as these algorithms improve, and co-located high-resolution cameras (such as the Multi-Angle Snowflake Camera instrument; Garrett et al., 2012) provide complementary information about selected individual particles, while the PIP, with its wider field of view, provides more information about the PSD and bulk snowfall properties.

Table 9 Mean bias, absolute error (MAE), and correlation coefficient (r) of PIP-derived versus D3R-estimated mean area-weighted ellipse ($\stackrel{\mathrm{‾}}{a_{\mathrm{ell}}}$) and box ($\stackrel{\mathrm{‾}}{a_{\mathrm{box}}}$) aspect ratio for each of the four retrieval experiments.

6 Conclusions

This study describes an algorithm that makes use of both polarimetric and dual-frequency radar measurements to retrieve microphysical properties of falling snow, including snowfall rate (volumetric and water equivalent); ice water content; particle size distribution; the relative contribution of pristine, aggregate, and rimed species; and particle orientation distribution. The algorithm is flexible in that it can use as many or as few measurements as available. In this study, it is applied to the Ku- and Ka-band measurements of the NASA D3R radar obtained during the ICE-POP 2018 field campaign but can be applied to additional or different frequencies. This is possible because it makes use of the ARTS microwave single-scattering property database for oriented particles (Brath et al., 2020), which encompasses ADDA (Yurkin and Hoekstra, 2011) scattering calculations over a wide range of frequencies. This differentiates it from methods that use T-matrix or Rayleigh–Gans approximations but constrains it to use the particle geometries that are available (at this time, only hexagonal plates and aggregates composed of these plates). More geometries are available for randomly oriented particles, but these cannot make use of the polarimetric information (although they are used in this study to represent the rimed particles).

The retrieval uses optimal estimation to solve for several parameters that describe the PSD, relative contribution of each species, and the orientation distribution along an entire radial simultaneously. This is necessary (versus a gate-by-gate approach) to account for the measurements sensitive to propagation effects (e.g., DWR and ϕ_dp). Examination of Jacobian matrices and averaging kernels shows that the DWR provides information regarding the characteristic particle size, and to a lesser extent, the rime fraction and shape parameter of the size distribution. The Z_dr measurements provide information regarding the mass fraction of pristine particles and their characteristic size and orientation distribution. Meanwhile, the ϕ_dp measurements are sensitive to most of the same measurements as Z_dr but are also sensitive to the overall particle concentration. Thus, by combining the dual-frequency and polarimetric measurements, some ambiguities can be resolved that should allow a better determination of the particle size distribution parameters and integrated quantities (e.g., ice water content, snowfall rate) than can be retrieved from single-frequency polarimetric measurements or dual-frequency, single-polarization measurements.

The D3R ICE-POP retrievals were validated using PIP and Pluvio measurements taken nearby at the May Hills ground site. The PIP measures the snow PSD directly (Pettersen et al., 2020) and several useful parameters can be derived directly from its measurements or indirectly with additional assumptions. These include the snowfall rate (volumetric and water equivalent), mean volume-weighted particle size, and effective density (Tokay et al., 2022), as well as parameters describing the mean aspect ratio and orientation distribution. We validated the retrieval during three events representing both warm and cold snow regimes (Kim et al., 2021). These events were chosen based upon availability of both PIP and D3R data, significant accumulation at the MHS location, and absence of any mixed-phase precipitation which the algorithm does not account for. Four retrieval experiments were performed to evaluate the utility of different measurement combinations: Ku-only, DWR-only, Ku-pol, and All-obs. In terms of mean absolute error and correlation, the volumetric snowfall rate was best retrieved (r=0.95), followed closely by the snow water equivalent rate (r=0.92). The Ku-pol method had the highest correlation to these parameters, while the DWR-only and All-obs methods had the lowest magnitude bias. These methods that incorporated DWR also had the best correlation to particle size (r=0.74), although none of the methods retrieved density particularly well (r=0.46). The ability of the measurements to retrieve mean aspect ratio was also inconclusive, although the polarimetric methods (Ku-pol and All-obs) had reduced biases and MAE relative to the Ku-only and DWR-only methods. The significant biases in particle size and snowfall rate appeared to be related to biases in the measured DWR (positive on 28 February and negative on 7–8 March), emphasizing the need for accurate DWR measurements and frequent calibration (e.g., co-located measurements at a non-attenuating frequency such as S or C band). Notwithstanding these calibration biases, during the most well-behaved event (7–8 March), where the PIP-derived reflectivity was closest to the D3R measurement, the All-obs method provided the best snowfall accumulation and closely approximated the observed time series of snowfall rate and particle size.

The D3R is scheduled to be deployed in Storrs, Connecticut, USA, during the 2021–2022 and 2022–2023 winters as part of the NASA-sponsored Investigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms (IMPACTS) field experiment. A similar deployment setup is planned with nearby PIP measurements. Additionally, the airborne Ku- and Ka-band HIWRAP radar will be available to evaluate the D3R calibration, and other lessons learned (e.g., video monitoring to observe snow accumulation on the radome, siting to avoid ground clutter) will be applied. Ongoing improvements to the PIP processing algorithms, particularly regarding the estimation of particle aspect ratio, will also be advantageous to further refine the algorithm described in this work. Availability of scattering databases for oriented particles with different geometries will facilitate running these retrievals as an ensemble to provide more robust posterior distributions of the retrieved parameters. Finally, the methodology can be expanded to accommodate liquid and melting particles, although scattering databases for the latter, particularly with the polarimetric parameters from oriented melting particles, are not yet mature enough for this application.

Appendix A: Derivation of single-scattering properties

Faithfully simulating the observables from polarimetric radar requires considering the incident and scattered Stokes vectors; these vectors are related via (Adams and Bettenhausen, 2012)

$$\begin{array}{}\text{(A1)}& \left[\begin{array}{c}{I}_{\mathrm{s}}\\ Q_{\mathrm{s}}\\ U_{\mathrm{s}}\\ V_{\mathrm{s}}\end{array}\right]={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}\left[\begin{array}{cccc}{Z}_{\mathrm{11}}& Z_{\mathrm{12}}& Z_{\mathrm{13}}& Z_{\mathrm{14}}\\ Z_{\mathrm{21}}& Z_{\mathrm{22}}& Z_{\mathrm{23}}& Z_{\mathrm{24}}\\ Z_{\mathrm{31}}& Z_{\mathrm{32}}& Z_{\mathrm{33}}& Z_{\mathrm{34}}\\ Z_{\mathrm{41}}& Z_{\mathrm{42}}& Z_{\mathrm{43}}& Z_{\mathrm{44}}\end{array}\right]\left[\begin{array}{c}{I}_{\mathrm{i}}\\ Q_{\mathrm{i}}\\ U_{\mathrm{i}}\\ V_{\mathrm{i}},\end{array}\right]\end{array}$$

where I, Q, U, and V are the elements of the Stokes vector, r is the distance from the sensor to the particle, Z_lm are the elements of the scattering or phase matrix, and the i and s subscripts indicate incidence and scattering, respectively.

To generate tables of the polarimetric, single-scattering properties, we transform the Stokes matrix elements to single-scattering properties more commonly used in radar meteorology, such as backscatter cross section at horizontal and vertical polarizations (${\mathit{\sigma }}_{\mathrm{hh}}^{\mathrm{b}}$ and ${\mathit{\sigma }}_{\mathrm{vv}}^{\mathrm{b}}$, respectively). Additionally, we include the complex, co-polar conjugate product C_hv between the scattering amplitude matrix elements at horizontal and vertical polarization (S_hh and S_vv, respectively) that allow for the co-polar correlation coefficient ρ_hv to be calculated. These variables are defined in terms of the phase matrix elements as (Ekelund et al., 2020)

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{hh}}^{\mathrm{b}}=\mathrm{2}\mathit{\pi }\left(Z_{\mathrm{11}}+Z_{\mathrm{12}}+Z_{\mathrm{21}}+Z_{\mathrm{22}}\right)\text{(A3)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{vv}}^{\mathrm{b}}=\mathrm{2}\mathit{\pi }\left(Z_{\mathrm{11}}-Z_{\mathrm{12}}-Z_{\mathrm{21}}+Z_{\mathrm{22}}\right)\text{(A4)}& {\displaystyle }\mathit{Re}\left(C_{\mathrm{hv}}\right)={\displaystyle \frac{Z_{\mathrm{44}}+Z_{\mathrm{33}}}{\mathrm{2}}}\text{(A5)}& {\displaystyle }\mathit{Im}\left(C_{\mathrm{hv}}\right)={\displaystyle \frac{Z_{\mathrm{43}}-Z_{\mathrm{34}}}{\mathrm{2}}}.\end{array}$$

All the expressions above define the phase matrix elements as in the backscatter direction, or opposite the incident direction, and the phase matrix elements have units of mm².

Expressions are given below for the radar reflectivity factor

$$\begin{array}{}\text{(A6)}& Z_{\mathrm{h},\mathrm{v}}={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K_{\mathrm{w}}{|}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}N\left(D\right)p\left(\mathit{\beta }\right){\mathit{\sigma }}_{\mathrm{h},\mathrm{v}}(D,\mathit{\beta })\mathrm{d}\mathit{\beta }\mathrm{d}D\end{array}$$

in units of mm⁶ m⁻³ and co-polar correlation coefficient

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& {\mathit{\rho }}_{\mathrm{hv}}={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}|K_{\mathrm{w}}{|}^{\mathrm{2}}}}\\ & {\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\int }_{\mathrm{0}}^{\mathit{\pi }}N\left(D\right)p\left(\mathit{\beta }\right)\sqrt{\mathit{Re}\left(C_{\mathrm{hv}}\right(D,\mathit{\beta }){)}^{\mathrm{2}}+\mathit{Im}{\left(C_{\mathrm{hv}}(D,\mathit{\beta })\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\beta }\mathrm{d}D}{\sqrt{Z_{\mathrm{h}}{Z}_{\mathrm{v}}}}}.\end{array}\end{array}$$

The propagation variables, K_dp, and the extinction cross sections (${\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{e}}$ and ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{e}}$) can be calculated from the corresponding extinction matrix of the particle. The oriented, azimuthally random uniform particles we use from the ARTS database have only three independent extinction matrix elements: K₁₁, K₁₂, and K₃₄. The propagation variables are thus defined as (Ekelund et al., 2020)

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{e}}=K_{\mathrm{11}}+K_{\mathrm{12}}\text{(A9)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{e}}=K_{\mathrm{11}}-K_{\mathrm{12}}\text{(A10)}& {\displaystyle }{K}_{\mathrm{dp}}={\displaystyle \frac{\mathrm{180}\times {\mathrm{10}}^{-\mathrm{3}}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}N\left(D\right)p\left(\mathit{\beta }\right)K_{\mathrm{34}}(D,\mathit{\beta })\mathrm{d}\mathit{\beta }\mathrm{d}D.\end{array}$$

The extinction matrix elements have units of mm² and K_dp has units of degree km⁻¹.

Appendix B: Radar forward model

The D3R measurements that comprise the observation vector Y (Eq. 12) are forward modeled from the measured Ku-band vertically polarized reflectivity $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{m}}$¹ and the state vector X. This combination of frequency and polarization was chosen because it is least prone to error in the attenuation correction. To simulate the D3R measurements, we must obtain both the backscattering properties of the ice particles as well as the propagation scattering properties (i.e., the propagation phase shifts and attenuation at each polarization) from the LUTs defined in Sect. 3. The components of Y are defined by the following equations that include the propagation effects along a radial:

$$\begin{array}{}\text{(B1)}& {\displaystyle }\text{DWR}\left(r\right)=\mathrm{10}{\log }_{\mathrm{10}}{Z}_{\mathrm{v},\mathrm{Ku}}^{\mathrm{m}}\left(r\right)-\mathrm{10}{\log }_{\mathrm{10}}{Z}_{\mathrm{v},\mathrm{Ka}}^{\mathrm{m}}\left(r\right),\text{(B2)}& {\displaystyle }{Z}_{\mathrm{dr}}(r,f)=\mathrm{10}{\log }_{\mathrm{10}}{Z}_{\mathrm{h},\mathrm{f}}^{\mathrm{m}}\left(r\right)-\mathrm{10}{\log }_{\mathrm{10}}{Z}_{v,f}^m\left(r\right),\text{(B3)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{dp}}(r,f)=\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{r}}{\int }}{K}_{\mathrm{dp}}(r^{\prime },f)\mathrm{d}r{}^{\prime }+{\mathit{\varphi }}_{\mathrm{sys}},\end{array}$$

where K_dp(f) is the specific differential phase at frequency f, ϕ_sys is the system differential phase upon transmission, and $Z_{\mathrm{p},\mathrm{f}}^{\mathrm{m}}$ is the attenuated reflectivity at polarization p and frequency f defined as

$$\begin{array}{}\text{(B4)}& Z_{\mathrm{p},\mathrm{f}}^{\mathrm{m}}\left(r\right)=\exp \left[-\mathrm{2}\underset{\mathrm{0}}{\overset{r}{\int }}{k}_{\mathrm{ext}}(p,f,r{}^{\prime })\mathrm{d}r{}^{\prime }\right]Z_{\mathrm{p},\mathrm{f}}^{\mathrm{e}}\left(r\right).\end{array}$$

Here, $Z_{\mathrm{p},\mathrm{f}}^{\mathrm{e}}$ is the intrinsic effective reflectivity factor and k_ext(p,f) is the specific attenuation, which includes contributions from gases, cloud liquid water, and hydrometeors. The symbols p and f indicate the polarization and frequency, respectively.

DWR, Z_dr, and ϕ_dp are forward simulated from Eqs. (B1)–(B4) by integrating numerically along each radial, following a conceptually similar process to that of Grecu et al. (2011). First, at each range gate, the attenuation contributions from gases (determined by the closest-in-time sounding from DGW) and cloud liquid (part of the retrieval state vector) are calculated. Next, $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{e}}$ is calculated by inserting the observed $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{m}}$ and the calculated two-way path-integrated attenuation from previous range gates into Eq. (B4). Then, $Z_{\mathrm{p},\mathrm{f}}^{\mathrm{e}}$, k_ext(p,f), and K_dp(f), along with the simulated PIP measurements in Table 4, are calculated from interpolated LUTs that have been combined for all the species according to the retrieval parameters and scaled to the observed $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{m}}$. These quantities are then used to forward model the measurements that comprise Y following Eqs. (B1)–(B4).

The process of simulating the radar measurements from combined LUTs can be described by considering a LUT for a variable V selected from Table 4. First, the LUT dimensionality for each species is reduced by linear interpolation over the following parameters: the radial zenith angle θ, temperature at the altitude of the range gate, and species-specific μ and κ values interpolated from the nodes for these parameters. At this stage, the remaining LUT dimensions are D_m and frequency (for the radar variables). Next, the unrimed aggregate (V_u) and rimed (V_r) LUTs are combined into a merged “aggregate” LUT (V_a) by weighting each species according to the rime mass fraction interpolated to each range gate r:

$$\begin{array}{}\text{(B5)}& V_{\mathrm{a}}\left(D_{\mathrm{m}}\right)=\left(\mathrm{1}-f_{\mathrm{r}}\left(r\right)\right)V_{\mathrm{u}}\left(D_{\mathrm{m}}\right)+f_{\mathrm{r}}\left(r\right)V_{\mathrm{r}}\left(D_{\mathrm{m}}\right).\end{array}$$

Recall from Sect. 3 that the LUTs are normalized to contain a constant 1 g m⁻³ ice water content; there exists a corresponding $N_{\mathrm{P}\mathrm{0}}^{*}\left(D_{\mathrm{m}}\right)$ value consistent with that water content. Because $N_{\mathrm{P}}^{*}$ is prescribed at each range gate from the interpolated retrieval parameter $\mathit{\delta }{N}_{\mathrm{P}}^{*}\left(r\right)$ (see definition in Table 5), the LUT for each species s is rescaled element wise to be consistent with this prescribed $N_P^{*}$:

$$\begin{array}{}\text{(B6)}& V_s^x\left(D_{\mathrm{m}}\right)=V_s\left(D_{\mathrm{m}}\right){\displaystyle \frac{N_{\mathrm{P}}^{*}\left(r\right)}{N_{\mathrm{0}}^{*}\left(D_{\mathrm{m}}\right)}}.\end{array}$$

The pristine and aggregate LUTs are then merged – individually for the scaled and unscaled LUTs – according to the interpolated retrieval parameters f_p (pristine mass fraction) and D_p (ratio of the pristine to aggregate D_m). To perform this merging, first the mean mass-weighted diameter dimension of the combined LUT (D_m,c) is calculated from the prescribed f_p, D_p, and aggregate D_m,a:

$$\begin{array}{}\text{(B7)}& D_{\mathrm{m},\mathrm{c}}=f_{\mathrm{p}}\left(r\right)D_{\mathrm{p}}\left(r\right)D_{\mathrm{m},\mathrm{a}}+\left(\mathrm{1}-f_{\mathrm{p}}\left(r\right)\right)D_{\mathrm{m},\mathrm{a}}.\end{array}$$

Next, the variables in the pristine LUTs (V_p and $V_{\mathrm{p}}^x$) are interpolated to this new set of D_m,c values and merged with the aggregate LUTs (V_a and $V_{\mathrm{a}}^x$) according to the prescribed mass fraction to create a combined LUT:

$$\begin{array}{}\text{(B8)}& {\displaystyle }{V}_{\mathrm{c}}\left(D_{\mathrm{m},\mathrm{c}}\right)=f_{\mathrm{p}}\left(r\right)V_{\mathrm{p}}\left(D_{\mathrm{p}}\right(r\left)D_{\mathrm{m},\mathrm{a}}\right)+(\mathrm{1}-f_{\mathrm{p}}(r\left)\right)V_{\mathrm{a}}\left(D_{\mathrm{m},\mathrm{a}}\right),\text{(B9)}& {\displaystyle }{V}_{\mathrm{c}}^x\left(D_{\mathrm{m},\mathrm{c}}\right)=f_{\mathrm{p}}\left(r\right)V_{\mathrm{p}}^x\left(D_{\mathrm{p}}\right(r\left)D_{\mathrm{m},\mathrm{a}}\right)+(\mathrm{1}-f_{\mathrm{p}}(r\left)\right)V_{\mathrm{a}}^x\left(D_{\mathrm{m},\mathrm{a}}\right).\end{array}$$

Now there exists a one-dimensional LUT for each parameter (and frequency for the radar variables) that is indexed by D_m and is consistent with the prescribed PSD characteristics from the retrieval state vector. The retrieved value of D_m(r) is that which produces the attenuation-corrected (from Eq. B4 over prior range gates) $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{c}}\left(r\right)$ from the $N_{\mathrm{P}}^{*}$-scaled Z_v,Ku LUT. From this retrieved D_m(r), the ice water content W(r) (in g m⁻³) is simply the scaling factor that is required to produce this $Z_{\mathrm{v},\mathrm{Ku}}^{\mathrm{c}}\left(r\right)$ from the unscaled LUT. Once D_m(r) and W(r) are determined, all of the remaining LUT parameters, including those needed to simulate the measurements in Y as well as the simulated PIP measurements used for validation, can be readily calculated by interpolating the LUTs to D_m(r) and scaling by W(r). These parameters are saved to the output data structure, and the propagation equations are iteratively integrated to repeat the process at the next range gate.