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 Kuand Ka-band measurements of the NASA Dual-polarization, Dual-frequency Doppler Radar (D3R) obtained during the International Collaborative Experiment for PyeongChang Olympic and Paralympics (ICE-POP 2018) field campaign, and incorporates the Atmospheric Radiative Transfer Simulator (ARTS) microwave single 5 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 show 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 10 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 mea15 surements 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 20 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 MAE relative to the 1 https://doi.org/10.5194/amt-2021-264 Preprint. Discussion started: 1 September 2021 c © Author(s) 2021. CC BY 4.0 License.

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 show 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 10 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 mea-15 surements 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 20 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 MAE relative to the 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 60 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, 65 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 70 with these capabilities deployed in locations subject to frequent and high-accumulation snowfall events. The NASA Dualpolarization, 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 75 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)  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 85 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 90 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 section 2, the particle scattering properties we use and the construction of lookup tables from these databases is described in section 3, algorithm mechanics and information content are analyzed in section 4, validation for selected cases is presented in section 5, and the conclusions are given in section 6. 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 100 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 789m msl). The D3R was configured to measure 150m range gates out to a maximum range of 39.75 km, combining a short pulse for ranges < 3.3km with a medium pulse for the remaining range gates. This gives a Ku-band 105 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 PPI scan at 5 • elevation follow by RHI scans at 51 • , 231 • , and 330 • azimuths, with each set of scans taking 5 minutes. 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 the D3R data available from the NASA ICE-POP data archive held at the Global Hydrology Resource Center 110 . 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. 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. The D3R Z dr and φ dp were calibrated on an event basis, and the absolution calibration of Ku and Ka Z was 115 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. 120 We examined time series of the pdfs of these quantities for the events listed in Table 1. In Figure 1, the top two rows show the time series of the pdf of Ku and Ka Z dr at elevation angles > 80 • and altitudes above 1.5km msl. 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 > 80 • is equal to zero. The time series of the calibrated Ku and Ka Z dr are shown in the third and fourth rows of Figure 1. In addition to the modification to Z dr , we also noticed some 125 unusual behavior in the small-particle DWR pdf (fifth row of Figure 1) during the February 28 case. When Z Ku < 10dBZ, -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 0600 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 (V. Chandrasekahr, pers. comm.), but for this study we have calibration from a ground radar.

PIP
The primary source of validation for the microphysical retrievals in this study is the Precipitation Imaging Package (PIP; Pettersen et al., 2020). 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 135 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 by 480 pixels with a resolution of 0.1 mm by 0.1 mm. The field of view (FOV), including the edge effects, is then 64 − D eq mm by 48 − D 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 117/D 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 140 particles with uniform size of D eq = 1 mm, the sampling volume is 790 m 3 for a one minute 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 145 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. This strategy uses the method implemented in the fit_ellipse program in the 150 Coyote IDL Library (http://www.idlcoyote.com) 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 minute 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 155 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.

160
Radiosonde launches were performed every three hours during precipitation events from the DGW Regional Weather Office, supplementing the normal 12-hourly observations, and used Meteomodem M10 radiosondes (Meteomodem, cited 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 supersatu-165 rated with respect to ice).

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 particle size distribution (PSD).

170
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) 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 GHz and 35.6 GHz), temperatures (190K, 230K, and 270K), 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 175 orientation angles (we refer to this angle hereafter as β). Table 2 lists values of these properties that correspond to the scattering calculations.
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 180 over a particle size distribution (PSD) using a modified gamma form (e.g., Petty and Huang, 2011): where D is the diameter of an equivalent-volume solid ice sphere in mm, D m is the mass-weighted mean equi-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 185 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 0 (D m , µ) is a normalized concentration factor such that the ice Table 3. Dimensions of the lookup tables derived from the particle types selected from the ARTS Single Scattering Database.
Dimension Range Orientation dispersion parameter (κ) 0-25 water content of all PSDs stored in the LUT is 1 g m −3 : where ρ i is the density of solid ice (0.917 g cm −3 ). The distribution is later scaled by a retrieved concentration factor to match 190 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: where κ is a dispersion parameter and G(κ) is a normalization factor such that the sum of probabilities is equal to one. When κ = 0, the distribution is uniform; as κ increases, the distribution becomes narrower and can be approximated by a normal distribution with standard deviation 1/κ. 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 200 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 Section 4.
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 205 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 copolar conjugate product of scattering Table 4. Variables stored within the lookup tables.

Group Variables Dimensions
Radar Variables Zvv, Z hh , kext,v, k ext,h , K dp , f,T,θ, Dm, µ, κ The formulas used to forward-model the radar measurements from these quantities are given in Appendix B.

210
In addition to the radar variables, we simulate several variables (denoted with subscript P ) that represent PIP-measured quantities. Because the PIP measures 2D projections of 3D 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 2D 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 215 to derive the 3D properties from the 2D 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 * P , area-weighted mean particle volume diameter D P , mean box aspect ratio (a box ) and mean ellipse aspect ratio (a ell ). The final two quantities are weighted by the projected area. The snowfall 220 water equivalent rate is calculated as: where m(D) is the mass of the particle in grams, V t (D, β) is the terminal velocity in m s −1 , and S is in units of mm hr −1 .
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 density/viscosity of air typical at the MHS 225 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: 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): The normalized intercept N * P is adapted from the definition given for N * 0,23 in Field et al. (2005): where d max2D (D, β) is the average 2D projected maximum particle dimension. This parameter is useful for providing temperaturedependent constraints on the PSD as will be shown in section 4.

235
D P is the ratio of the 4th to 3rd moments of the PSD in terms of D eq : 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 (a ell ): where b a (D, β) 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 (a box ): where y x (D, β) is the ratio of the vertical to horizontal dimensions of the bounding box of the simulated PIP image.

245
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 Figure 2, 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 < 1mm and D P < 3mm, 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 250 been observed coincident with large aggregates (20mm projected diameter) in dendritic growth regimes during OLYMPEX , 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 255 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 κ ( Figure 3). 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 260 plates also does not depend much on size, although some resonance effects lead to a small decrease in K dp at larger size. From Figures 2 and 3, 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.

Algorithm Description
Optimal estimation (OE; Rodgers, 2000) is a form of Bayesian inversion that assumes Gaussian error statistics and accommo-265 dates 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.
270 Figure 3. Polarimetric variables Z dr (top) and K dp (bottom) as a function of mean particle size Dm and orientation distribution parameter κ. The left column represents aggregates, and the right column is for individual plates. All values are for horizontal incidence at Ku band at 270 K.

Optimal estimation setup
The cost function that is minimized by optimal estimation is 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 275 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: 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 280 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:

285
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 600m horizontally and 250m 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.

290
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 under-utilization 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 295 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 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 accomodate lognormal distributions and/or to linearize the forward model.   (Schrom et al., 2021). Although these forward model errors are correlated between variables, analysis of the covariance 300 matrices showed the correlation between DW R 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.
The optimal estimation process along a ray of radar data is illustrated in Figure 4, which shows the observed (Y obs ) and simulated (Y sim ) measurements, and Figure 5, 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 a DWR peaks at 0-5 km 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 km and 30 km range, bracketing the DWR peak. The φ dp increases most rapidly in the 10-20 km range with 310 smaller rates of increase elsewhere, reaching 15 • at Ku band and 50 • at Ka band. The initial profiles of the retrieval parameters N * 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 Figure 5. The iterative adjustments guided by the Jacobian matrix respond to the initial differences between Y obs and Y sim : -N * P decreases slightly near the DWR peaks and increases elsewhere;

315
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;

320
µ 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 imply that there is little differential attenuation, and the observed near-field DWR can be attributed to particle size effects.

325
-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 .
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 330 histogram for this case in Figure 1 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.

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 is 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:

340
-Increasing N * 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 contribute more to 345 reflectivity, but decreases φ dp downrange as the ice water content in 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 no 350 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 355 a requirement for reducing ambiguity, or cross-talk, in the retrieved state. The information content and cross-talk among Figure 6. 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 Sy -i.e., the expected standard deviation of that measurement's error. 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: Values close to one indicate strong influence of measurements, and values close to zero indicate that the retrieval is heavily 360 influenced by the prior. In Figure 7, 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 * P ), which both the DWR and φ dp are highly sensitive to ( Figure 6). 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 * P , owing to their sensitivity to the polarimetric parameters, but a lower 10% quantile, likely originating from situation where the f r is 365 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% quantile 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 370 significant cross-talk, i.e., a strong correlation in the posterior state vector, between the following groups of variables: (N * P , Figure 7. Quantiles of the averaging kernel matrix diagonal elements and off-diagonal elements representing different parameters at the same range gate.
µ 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 ( Figure 3, and the shape parameter (µ) of the pristine PSD is not going to influence 375 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 one. 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 380 values will induce large changes to the DWR. However, there were several RHI scans, particularly on February 28, where some cloud liquid was necessary to explain high DWR values that could not be achieved by differential scattering alone.

Validation
The primary tool used for validation of the D3R retrievals is the PIP located at the MHS site. 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 385 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: Frequency Precipitation Radar (DPR);

390
-Ku-pol: A single-frequency polarimetric retrieval using Z dr and φ dp at 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., 2021). Although there will be some departures from Mie theory for non-spherical particles, at Ku band these are relatively small (Kuo et al., 2016). The larger 395 contribution to error is the various assumptions required to derive particle density from the PIP size and fall speed (Tokay et al., 2021). 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 Figure 8).
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 MRR 400 radar profiles indicate significant echo growth (perhaps due to aggregation) at low levels, which may have continued in the clutter region. The 28 Feb case exhibits a similar bias between the PIP-derived and D3R-measured reflectivity at Ku band in the first 6 hours, after which the D3R reflectivity comes within the lower range of PIP estimates. This was a deeper, warm low case, and the MRR and 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 0400 UTC and the air 405 temperature was above 0 • C until 0800 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 section 2). The 7-8 Mar 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 • C 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 410 listed in Table 4. Because the PIP measurements are taken every minute, whereas the D3R RHI scans were conducted every 5-6 minutes, 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 minutes, depending on 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 415 to ensure that the PIP measurements are reasonably representative of the D3R observations.

Snowfall rate and water equivalent
The time series of snow water equivalent rate (S) for each event are shown in Figure 9. The D3R substantially underestimates snowfall throughout the 9 Jan and 28 Feb events, which is not surprising since the PIP-derived reflectivity significantly exceeds Z dr (thrid 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., 2021), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.
the D3R measurement for reasons discussed previously in this section. It is worth noting, however, that the DWR-only and 420 All-Obs retrievals are in good agreement with the Pluvio-measured S accumulation on 9 Jan, which would be consistent with aggregation processes in the lowest levels that increase reflectivity, but do not increase S. The 7-8 Mar 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 425 DWR time series for this event in Figure 8 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 Jan event higher than the retrievals and in better agreement with the PIP. calculations involve a variety of assumptions to derive particle mass from shape and fall speed (Tokay et al., 2021), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.
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 Feb and 7-8 Mar events, and the MAE is highest for this method. The Ku-pol method provides the best correlation for S V and S even though it has the 435 largest magnitude bias of all the experiments.

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 Figure 10. On 9 Jan 2018, all of the radar retrievals were biased low 440 with respect to the PIP, which is consistent with the reflectivity bias we noted on this day. On 28 Feb, the PIP measured a rapid increase of D P to a maximum at 06 UTC, followed by a decrease to a minimum around 09 UTC and another maximum at 15 UTC. None of the retrievals did a particularly good job of capturing the peak at 06 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 Mar, Table 7. Mean bias, absolute error (MAE), and correlation coefficient (r) of PIP-derived versus D3R-estimated volumetric snowfall rate (SV ) 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.  Figure 10. Time series of D3R-and PIP-derived mean volume-weighted particle diameter (top) and effective density (bottom) for each event.
Because the PIP calculations involve a variety of assumptions to derive particle density from shape and fall speed (Tokay et al., 2021), a range between the minimum and maximum from these calculations is shaded around the mean PIP-derived value.
both methods that used the DWR (DWR-only and All-obs) captured the temporal variability of D P quite well, but were biased 445 low; this is consistent with a suspected low bias in the DWR for this event (mean cloud-top DWRs were slightly below zero; see Figure 1).
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 counter-intuitive, since this retrieval does not consider 450 the DWR, which is the measurement most sensitive to D P . However, the suspected low bias of the DWR on 7-8 Mar, 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 455 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. 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 Figure 10 shows the time series of ρ P for each event. In general, we find that when 460 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 Feb, all of the D3R methods are in tight agreement whereas on 7-8 Mar, 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 465 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).

Bulk particle orientation and aspect ratio 470
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 (a ell ) and box (a box ) aspect ratios. In Figure 11, the influence of the Z dr and K dp measurements (Figure 8) 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 UTC on 9 Jan result in low aspect ratios at that time.

475
There is not much of a trend in Z dr on 28 Feb, owing to the constant presence of large aggregates which dominate the Z dr measurement, but there is a notable peak in K dp around 09 UTC which corresponds to the minimum aspect ratio for this event.
The 7-8 Mar event had the least variable polarimetric signatures, but a short peak in Z dr around 18 UTC on 7 Mar corresponds to a drop in aspect ratio at that time.
The comparison of the retrieved aspect ratios to those derived from the PIP is inconsistent from event to event. There appears 480 to be little variability in the PIP time series on 9 Jan for either measure of aspect ratio. There is a steady increase in a ell between 09 and 12 UTC on 28 Feb, matching the retrieved behavior, although the PIP range is considerably smaller than the retrieved range. The 7-8 Mar 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 a ell ) and the DWR-only 485 method (for a box ). This is a surprising result, especially in light of the very small range of these retrieved values in Figure 11.
In these methods, the very limited AR 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.

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 495 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 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 500 geometries are available for randomly-oriented particles, but these can not 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., DW R and φ dp ). Examination of Jacobian ma-505 trices and averaging kernels show 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 510 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), these vectors are related via (Adams and Bettenhausen, 2012) 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-555 scattering properties more commonly used in radar meteorology, such as backscatter cross-section at horizontal and vertical polarizations (σ b hh and σ b vv , 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 copolar correlation coefficient ρ hv to be calculated. These variables are defined in terms of the phase matrix elements as 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 2 .
Expressions are given below for the radar reflectivity factor The propagation variables, K dp and the extinction cross-sections (σ e h and σ e v ) 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 11 , K 12 , and K 34 . The propagation variables are thus defined as  575 The extinction matrix elements have units of mm s and K dp has units of • km −1 .

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 m v,Ku 1 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 585 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 lookup tables (LUTs) defined in Section 3. The components of Y are defined by the following equations that include the propagation effects along a radial: Z dr (r, f ) = 10log 10 Z m h,f (r) − 10log 10 Z m v,f (r), (B2) 590 φ dp (r, f ) = 2 r 0 K dp (r , f )dr + φ sys , from gases, cloud liquid water, and hydrometeors. The symbols p, f indicate the polarization and frequency, respectively.
DWR, Z dr , and φ dp are forward simulated from equations 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 e v,Ku is calculated by inserting the observed Z m v,Ku and the calculated 2-way path integrated attenuation from previous 600 range gates into Eq. B4. Then, Z e p,f , 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 m v,Ku . These quantities are then used to forward-model the measurements that comprise Y following equations B1-B4.
The process of simulating the radar measurements from combined LUTs can be described by considering a LUT for a 605 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: 610 V a (D m ) = (1 − f r (r))V u (D m ) + f r (r)V r (D m ).
Recall from Section 3 that the LUTs are normalized to contain a constant 1 g m −3 ice water content, there exists a corresponding N * P 0 (D m ) value consistent with that water content. Because N * P is prescribed at each range gate from the interpolated retrieval parameter δN * P (r) (see definition in Table 5), the LUT for each species s is re-scaled element-wise to be consistent with this prescribed N * P : The pristine and aggregate LUTs are then merged -individually for the scaled and unscaled LUTs -according to the interpo- 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 c v,Ku (r) from the N * P -scaled Z v,Ku LUT. From this retrieved D m (r), the ice water content W (r) (in g m −3 ) is simply the scaling factor that is required to produce this 630 Z c v,Ku (r) 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.
Author contributions. Dr. Munchak led this research by developing the lookup tables and optimal estimation methodology, processing the radar data, and calculating validation statistics. Dr. Schrom assisted with the generation of lookup tables, polarimetric components of the radar forward model, and assessment of forward model error. Dr. Tokay processed the PIP and Pluvio measurements that were used for validation. Dr. Helms performed additional processing of the PIP data to derive the bulk aspect ratio and orientation parameters.
Competing interests. S. J. Munchak is a member of the editorial board of Atmospheric Measurement Techniques. Otherwise, the authors declare that they have no conflict of interest.