Dual-Wavelength Radar Technique Development for Snow Rate Estimation: A Case Study from GCPEx

Quantitative Precipitation Estimation (QPE) of snowfall has generally been expressed in power-law form between equivalent radar reflectivity factor (Ze) and liquid equivalent snow rate (SR). It is known that there is large variability in the pre-factor of the power law due to changes in particle size distribution (PSD), density, and fall velocity whereas the variability of the exponent is considerably smaller. The dual-wavelength radar reflectivity ratio (DWR) technique can improve SR accuracy by estimating one of the PSD parameters (characteristic diameter) thus reducing the variability due to 15 the pre-factor. The two frequencies commonly used in dual-wavelength techniques are Kuand Ka-bands. The basic idea of DWR is that the snow particle size-to-wavelength ratio is such as to fall in the Rayleigh region at Ku-band but in the Mie region at Ka-band. We propose a method for snow rate estimation by using NASA D3R radar DWR and Ka-band reflectivity observations 20 collected during a long-duration synoptic snow event on 30–31 January 2012 during the GCPEx (GPM Cold Season Precipitation Experiment). Since the particle mass can be estimated using 2D-video disdrometer (2DVD) fall speed data and hydrodynamic theory, we simulate the DWR and compare directly with D3R radar measurements. We also use the 2DVDbased mass to compute the 2DVD-based SR. Using three different mass estimation methods, we arrive at three respective sets of Z-SR and SR(Zh, DWR) relationships. We then use these relationships with D3R measurements to compute radar25 based SR. Finally, we validate our method by comparing the D3R radar-retrieved SR with accumulated SR directly measured by a well-shielded Pluvio gauge for the entire synoptic event.


Introduction
A detailed understanding of the geometric, microphysical, and scattering properties of ice hydrometeors is a vital prerequisite for the development of radar-based quantitative precipitation estimation (QPE) algorithms.Recent advances in surface and airborne optical imaging instruments and the wide proliferation of dual-polarization and multi-wavelength radar systems (ground based, airborne or satellite) have allowed for observations of the complexity inherent in winter precipitation via dedicated field programs (e.g., Skofronick-Jackson et al., 2015;Petäjä et al., 2016).These large field programs are vital given that the retrieval problem is severely underconstrained due to large number of geometrical and microphysical parameters of natural snowfall, their extreme sensitivity to subtle changes in environmental conditions, and co-existence of different populations of particle types within the sample volume (e.g., Szyrmer and Zawadzki, 2014).
The surface imaging instruments that give complementary measurements and are used in a number of recent studies include (i) 2-D video disdrometer (2DVD; Schönhuber et al., 2008), (ii) precipitation imaging package (PIP; von Lerber et al., 2017), (iii) Multi-Angle Snowflake Camera G.-J. Huang et al.: Dual-wavelength radar technique development for snow rate estimation (MASC; Garrett et al., 2012).When these instruments are used in conjunction with a well-shielded GEONOR or PLU-VIO gauge, it is shown that a physically consistent representation of the geometric, microphysical, and scattering properties needed for radar-based QPE can be achieved (Szyrmer and Zawadzki, 2010;Huang et al., 2015;von Lerber et al., 2017;Bukovčić et al., 2018).In this study, we use the 2DVD and PLUVIO gauge located within a double fence international reference (DFIR) wind shield to reduce wind effects.
Radar-based QPE has generally been based on Z e -SR (Z e is reflectivity; SR is liquid equivalent snow rate) power laws of the form Z e = α(SR) β , where the prefactor and exponent are estimated based on (i) direct correlation of radarmeasured Z e with snow gauges (Rasmussen et al., 2003;Fujiyoshi et al., 1990;Wolfe and Snider, 2012) or (ii) using imaging disdrometers such as 2DVD or PIP (Huang et al., 2015;von Lerber et al., 2017).Recently, Falconi et al. (2018) developed Z e -SR power laws at three frequencies (X-, Ka-, and W-band) by direct correlation of radar and PIP observations.These studies have highlighted the large variability of α due to particle size distribution (PSD), density, fall velocity, and dominant snow type, whereas the variability in β is considerably smaller.Similarly, both methods, (i) and (ii), have been used to estimate ice water content (IWC) from Z e using power laws of the form Z e = a(IWC) b based on airborne particle probe data, direct measurements of IWC, and airborne measurements of Z e (principally at X-, Ka-, and W-bands) (e.g., Heymsfield et al., 2005Heymsfield et al., , 2016;;Hogan et al., 2006).The advantage of airborne data is that a wide variety of temperatures and cloud types can be sampled (Heymsfield et al., 2016).
The dual-wavelength reflectivity ratio (DWR, the ratio of reflectivity from two different bands) radar-based QPE was proposed by Matrosov (1998), Matrosov et al. (2005) to improve SR accuracy by estimating the PSD parameter (median volume diameter D 0 ) with relatively low dependence on density if assumed constant.There has been limited use of dual-λ techniques for snowfall estimation, mainly using vertical-pointing ground radars or nadir-pointing airborne radars (Liao et al., 2005(Liao et al., , 2008(Liao et al., , 2016;;Szyrmer and Zawadzki, 2014;Falconi et al., 2018).The dual-λ method is of interest to us due to the availability of the NASA D3R scanning radar (Vega et al., 2014), which, to the best of our knowledge, has not been exploited for snow QPE to date.
The DWR is defined as the ratio of the equivalent radar reflectivity factors at two different frequency bands.The main principle in DWR is that the particle's size-to-wavelength ratio falls in the Rayleigh region at a low-frequency band (e.g., Ku-band) but in the Mie region at a high-frequency band (e.g., Ka-band) (Matrosov, 1998;Matrosov et al., 2005;Liao et al., 2016).Previous studies have shown that the DWR can be used to estimate D m , where D m is defined as the ratio of the fourth moment to the third moment of the PSD expressed in terms of liquid-equivalent size or mass (Liao et al., 2016).In this sense the DWR is similar to differential reflectivity (Z dr ) in dual-polarization radar technique, where Z dr is used to estimate D m (but the physical principles are, of course, different; Meneghini and Liao, 2007).The SR is obtained by "adjusting" the coefficient α in the Z e -SR power law based on the estimation of D m provided by the DWR.The prefactor α depends on the intercept parameter of the PSD (von Lerber et al., 2017) and not on D m directly.However, because of the apparent negative correlation between D m and PSD intercept parameter for a snowfall of a given intensity (Delanoë et al., 2005;Tiira et al., 2016), measurements of D m can be used to adjust the Z e -SR power law.
This paper is organized as follows.In Sect.2, we introduce the approach and methodologies proposed and used in this study, which may be considered technique development.We briefly explain how to estimate the mass of ice particles using a set of aerodynamic equations based on Böhm (1989) and Heymsfield and Westbrook (2010).We also give a brief introduction of the scattering model based on particle mass.Section 3 provides a brief overview of instruments installed at the test site and the dual-wavelength radar used in this study (D3R: Vega et al., 2014).We analyze surface and D3R radar data from one synoptic snowfall event during GCPEx and compare SR retrieved from DWR-based relations with SR measured by a snow gauge.The conclusions and possibilities for further improvement of the proposed techniques are discussed in Sect. 4. The acronyms and symbols are listed in Appendix.

Estimation of particle mass
The direct estimation of the mass of an ice particle is difficult and at present there is no instrument available to do this automatically.The conventional method is to use a power-law relation between the mass and the maximum dimension of the particle of the form m = aD b , where the prefactor a and exponent b are computed via measurements of particle size distribution N (D) from aircraft probes and independent measurements of the total ice water content as an integral constraint (Heymsfield and Westbrook, 2010).A similar method was used by Brandes et al. (2007), who used 2DVD data for N (D) and a snow gage for the liquid equivalent snow accumulation over periods of 5 min.These methods are more representative of an average relation when one particle type (e.g., snow aggregates) dominates the snowfall with large deviations possible for individual events with differing particle types (e.g., graupel).
To overcome these difficulties a more general method was proposed by Böhm (1989) based on estimating mass from fall velocity measurements, geometry, and environmental data if the measured fall velocity is in fact the terminal velocity (i.e., in the absence of vertical air motion or turbulence and in more or less uniform precipitation).The methodol-ogy has been described in detail by Szyrmer and Zawadzki (2010), Huang et al. (2015), andvon Lerber et al. (2017), and we refer to these articles for details.The essential feature is the unique nonlinear relation between the Davies (1945) number (X) and the Reynolds number (Re), where X is the ratio of mass to area or m/A 0.25 r (A r = A e /A is the area ratio, where A e is the effective projected area normal to the flow and A is the area of the minimum circumscribing circle or ellipse that completely contains A e ) and the Re is the product of terminal fall speed and the characteristic dimension of the particle.We have neglected the environmental parameters (air density, viscosity) as well as boundary layer depth of Abraham (1970) and the inviscid drag coefficient.The procedure is to (i) compute Re from fall velocity measurements and characteristic dimension of the particle (usually the maximum dimension), (ii) compute the Davies number X, which is expressed as a nonlinear function of Re, and boundary layer parameters (C 0 = 0.6 and δ 0 = 5.83; Böhm, 1989) and (iii) estimate particle mass from X and A r .Heymsfield and Westbrook (2010) proposed a simple adjustment (based on field and tank experiments) by defining a modified Davies number as proportional to m/A 0.5 r along with different boundary layer constants (C 0 = 0.292; δ 0 = 9.06) from Böhm.Their adjustment was shown to be in very good agreement with recent tank experiments by Westbrook and Sephton (2017), especially for particles like pristine dendrites with low A r and at low Re.Note that the difference of C 0 and δ 0 in Böhm and Heymsfield-Westbrook equations is mainly due to differences in the shape-correcting factor (A r ) to find the optimal relation between drag coefficient (or Davies number, X) and Reynolds number (Re).This is the main parameterization error in this set of equations.

Geometric and fall speed measurements
One source of uncertainty in applying the Böhm or Heymsfield and Westbrook (HW) method is calculating the area ratio (A r ) using instruments such as 2DVD or precipitation instrument package (PIP) as they do not give the projected area normal to the flow (i.e., they do not give the needed top view, but rather the 2DVD gives two side views on orthogonal planes as illustrated in Fig. 1).This is reasonable for snow aggregates which are expected to be randomly oriented.The other source of uncertainty is in the definition of characteristic dimension used in Re, which in the HW method is taken to be the diameter of the circumscribing circle that completely encloses the projected area, the maximum dimension (D max ; this is what we use for the 2DVD in our application of the HW method).For the Böhm method we use the procedure in Huang et al. (2015), which used D app defined as the equalvolume spherical diameter.
The two-dimensional video disdrometer (2DVD) used herein is described in Schönhuber et al. (2000), and calibration and accuracy of the instrument are detailed in Bernauer et al. (2015).The 2DVD is equipped with two line-scan cam-eras (referred to as cameras A and B) which can capture the particle image projection in two orthogonal planes (two side views).As mentioned earlier the area ratio (A r ) should be obtained from the projected image in the plane normal to the flow (i.e., top or bottom view).However, to the best of our knowledge, there are no ground-based instruments that can automatically and continuously capture the horizontal projected views (i.e., in the plane orthogonal to the flow) of precipitation particles (however, 3-D-reconstruction based on multiple views can give this information; Kleinkort et al., 2017).Compared with other optical-based instruments, such as HVSD (hydrometeor velocity size detector; Barthazy et al., 2004) or SVI (snow video imager; Newman et al., 2009), which only captures the projected view in one plane, the 2DVD offers views in two orthogonal planes, giving more geometric information.Figure 1 shows a snowflake observed by a 2DVD from two cameras.The thick black line is the contour of the particle and the thin black lines show the holes inside the particle.The effective projected area A e in the definition of area ratio is easy to compute by counting total pixels from the particle's image, and then multiplying by horizontal and vertical pixel width.The blue line is the minimum circumscribed ellipse.The area of the ellipse is A in the definition of area ratio.The size of particle measured by 2DVD is called the apparent diameter (D app ) which is defined as the diameter of the equivalent volume sphere (Schönhuber et al., 2000;Huang et al., 2015).The D app is used when computing Re, as mentioned earlier.The area ratio and D app are the geometric parameters that are used in our implementation of the Böhm method.
In our application of the HW method, the A is based on the diameter of the circumscribed circle that completely encloses the projected pixel area (A e ), which is easy to calculate from the contours in Fig. 1.Thus the area ratio is A e /A, while the characteristic dimension in Re is the diameter of the circumscribing circle.Note that the area ratio and characteristic dimension in Re depend on the type of instrument used (e.g., advanced version of snow video imager by von Lerber et al., 2017;the HVSD by Szyrmer and Zawadzki, 2010).These instruments give a projected view in one plane only and thus geometric corrections are used as detailed in the two references.
The two optic planes of the 2DVD are separated by around 6 mm and the accurate distance is based on calibration by dropping 10 mm steel balls at three corners of the sensing area (details of the calibration as well as accuracy of size, fall speed, and other geometric measures are given in Bernauer et al., 2015).During certain time periods, more than one precipitation particle falls in the 2DVD observation area.Since the two cameras look in different directions, the particles observed by camera A and camera B need to be paired.This pairing procedure is called "matching", and it is illustrated in Fig. 2. The time period [t 1 , t 2 ] is dependent on the assumed reasonable fall speed range.Assuming that the minimum and maximum reasonable fall speeds are v min and v max , respec- tively, the distance between two optic planes is D d , and camera A observed a particle at t 0 , we have t 1 = t 0 + D d /v max and t 2 = t 0 + D d /v min .After matching, the fall speed can be calculated as D d / t, where t is the time difference between two cameras observing the same particle.Because the fall speed of the 2DVD is dependent on matching, the geometric features and fall speeds will be in error when mismatch occurs.Huang et al. (2010) analyzed snow data from the 2DVD and found that the 2DVD manufacturer's matching algorithm for snow resulted in a significant mismatching problem (see also Bernauer et al., 2015).In the Appendix of Huang et al. (2010), they showed that the mismatch will cause the volume, vertical dimension, and fall speed of particles to be overestimated.Subsequently, the mass of particles will also be overestimated, mainly because of fall speed.To get the best estimation of mass, they used 2DVD singlecamera data and re-did the matching based on a weighted Hanesch criteria (Hanesch, 1999).If the match criteria are not satisfied, then that particle is rejected; it follows that the concentration will tend to be underestimated.To readjust the measured concentration for this underestimate (assumed to be a constant factor), the procedure described in Huang et al. (2015) is used, which only involves the ratio of the total number of particles counted in the scan area of the single camera to the number of successfully matched particles in the virtual measurement area.For the event analyzed here (using method 1 in Sect.3.3), this adjustment factor is between 1.1 and 1.5.The Pluvio gauge accumulation is not used as a constraint in method 1.The disadvantage of using single-camera data, as described in Huang et al. (2015), is that the particle contour data are not available (i.e., the manufacturer's code does not provide line scan data from single camera).Without contour data, both D app and A can only be estimated by the maximum width of the scan line and height of the particle as detailed in Huang et al. (2015).Moreover, the diame-Figure 2. Illustration of the matching procedure.In the situation shown, it is assumed that camera A observed a particle at time t 0 , and afterwards during a certain time period, t 1 to t 2 , camera B observed two particles.The matching procedure decides which particle observed by camera B is the same particle observed by camera A.
ter of the circumscribing circle or ellipse cannot be obtained without contour data.The only quantity included in singlecamera data is A e in terms of number of pixels.The Huang and Bringi approach (Huang et al., 2015) is referred to as HB, because both PSD (particle size distribution) and reflectivity (Z e ) are computed using D app as the measure of particle size.
For methods 2 and 3 in Sect.3.3, we used the manufacturer's matching algorithm, which gives the contour data.To avoid overestimating mass due to mismatch, we need to filter out those particles with unreasonable fall speeds.The vertical dimension of the particle's image before match is expressed as a number of scan lines (i.e., how many scan lines are masked by the particle).After match (so v t is known), the vertical pixel width is v t /f s , where f s is the scan frequency of a camera (∼ 55 kHz), and the vertical size of the particle is the vertical pixel width multiplied by the number of scan lines.Because two optical planes of the 2DVD are parallel, theoretically, the number of scan lines from cameras A and B should be the same.Considering the distance from the particle to the two cameras (projective effect of a camera), the digital error of a camera and particle rotation in two planes, the difference in the number of scan lines between two cameras may not always be the same but should be very close.Hanesch (1999) gave a set of matching criteria, the most important being the tolerance of the number of scan lines between the two cameras (see Table 1).To obtain reliable fall speeds, we examined all matched particles (given by the manufacturer's matching algorithm, numbering 507 833 for the event and marked as green in Fig. 3) and removed those particles which did not satisfy the Hanesch scan line criteria, resulting in 175 199 (34.5 %) of particles that did satisfy the match criteria (magenta marks in Fig. 3).We used the fall speeds of these filtered particles to compute their mass for both Böhm and Heymsfield-Westbrook methods, and then divided the mass by apparent volume (= π D 3 app /6) to get the particle density.Since the maximum density of ice particles is around 0.9 g cm −3 , we further remove particles with density larger than 1 g cm −3 .After this two-step filtering, the particles we use for further analysis (numbering 128 063) are shown in Fig. 3 as blue points.The filtering will eliminate particles, which will reduce the liquid equivalent snow accumulation.Hence, the Pluvio gauge accumulation is used as an integral constraint, i.e., the concentration in each bin is increased by a constant factor to match the 2DVD accumulation to the Pluvio accumulation.This constraint is only used in methods 2 and 3 in Sect.3.3.

Scattering model
The scattering computation of ice particles is difficult because of their irregular shapes with large natural variability (e.g., snow aggregates or rimed crystals).The most common scattering method used in the meteorological community is the discrete dipole approximation (DDA; Draine and Flatau, 1994).However, DDA is very time consuming and not suitable for large numbers of particles, especially at Wband (e.g., Chobanyan et al., 2015).On the other hand, the T-matrix method (Mishchenko et al., 2002) is more time efficient and commonly used in radar meteorology but it requires that the irregular particle shape be simplified to an axis-symmetric shape (e.g., spheroid).Ryzhkov et al. (1998) have shown that, in the Rayleigh region, the radar cross section is mainly related to particle's mass squared and less to the shape.For Mie scattering, however, the irregular snow shape plays a more significant role.Westbrook et al. (2006Westbrook et al. ( , 2008) ) used the Rayleigh-Gans approximation to develop an analytical equation for the scattering cross sections of simulated snow aggregates of bullet rosettes using an empirical fit to the form factor that accounts for deviations from the Rayleigh limit.Here, we use two scattering models, one based on the soft spheroid (Huang et al., 2015) with a fixed axis ratio and quasi-random orientation.The apparent density is calculated as the ratio of mass to apparent volume.There is considerable controversy in the literature on the applicability of the soft spheroid model with a fixed axis ratio, especially at Ka and higher frequencies such as W-band (e.g., Petty and Huang, 2010;Botta et al., 2010;Leinonen et al., 2012;Kneifel et al., 2015).However, Falconi et al. (2018) used the soft spheroid scattering model using T-matrix to compute Z e (at X-, Ka-, and W-bands) and showed that an effective optimized axis ratio of (oblate) spheroid could be selected that directly matches measured Z e by radar (their optimal axis ratio, however, varied with the frequency band, i.e., 1 for X-band, 0.8 for Ka, and 0.6 for W).They also found some differences in the optimal axis ratios for fluffy snow vs. rimed snow.Nevertheless, they compared DDA calculations of complex-shaped aggregates to the soft spheroid www.atmos-meas-tech.net/12/1409/2019/Atmos.Meas.Tech., 12, 1409-1427, 2019 model at W-band and concluded that the axis ratio can be used as a tuning parameter.They also showed the importance of size integration to compute Z e , i.e., the product of N(D) and the radar cross section for the soft spheroid vs. complex-shape aggregates.Their result implied that smaller particles had a larger value for the product when using a soft spheroid of 0.6 axis ratio relative to complex aggregates and vice versa for larger particles, leading to compensation when Z e is computed by size integration over all sizes.Thus, the soft spheroid model with axis ratio at 0.8 used by Huang et al. (2015), and which is used herein at Ku-and Ka-bands, is a reasonable approximation.
The second scattering model we used herein is from Liao et al. (2013), who use an effective fixed density approach to justify the oblate spheroid model.To compare the scattering properties of a snow aggregate with its simplified equal-mass spheroid, Liao et al. (2013) used six-branch bullet rosette snow crystals with maximum dimensions of 200 and 400 µm as two basic elements that simulate snow aggregation.They computed the backscattering coefficient, extinction coefficient, and asymmetry factor for simulated snowflakes, using the DDA and for the corresponding spheres and spheroids with the same mass but density fixed at 0.2 or 0.3 g cm −3 , and hence the apparent sphere volume equals the mass divided by the assumed fixed density.They showed that, when the frequency was lower than 35 GHz (Ka-band), the Mie scattering properties of spheres with a fixed density equal to 0.2 g cm −3 were in a good agreement with the scattering results for the simulated complex-shaped aggregate model with the same mass using the DDA (see also Kuo et al., 2016).They also showed this agreement with a spheroid model with a fixed axis ratio of 0.6 and random orientation.Here, we use the Liao et al. (2013) equivalent spheroid model with a fixed effective density of 0.2 g cm −3 at Ku-and Ka-bands (note that we estimate the mass of each particle from 2DVD measurements as described in Sect.2.1).Note that this fixeddensity spheroid scattering model is not based on microphysics (where the density would fall off inversely with increasing size) but on scattering equivalence with a simulated (same-mass) complex-shaped aggregate snowflake (Liao et al., 2016).
3 Case analysis

Test site instrumentation and the synoptic event
The GPM Cold-season Precipitation Experiment (GCPEx) was conducted by the National Aeronautics and Space Administration (NASA), USA, in cooperation with Environment Canada in Ontario, Canada from 17 January to 29 February 2012.The goal of GCPEx was ". . . to characterize the ability of multi-frequency active and passive microwave sensors to detect and estimate falling snow. . ." (Skofronick-Jackson et al., 2015).The field experiment We examine a snowfall event on 30-31 January 2012 that occurred across the GCPEx study area between roughly 22:00 UTC, 30 January and 04:00 UTC, 31 January.Details of this case using King City radar and aircraft spiral descent over the CARE site is given in Skofronick-Jackson et al. (2015).This event resulted in liquid accumulations of roughly 1-4 mm across the GCPEx domain with fairly uniform snowfall rates throughout the event.At the CARE site the accumulations over an 8 h period were < 3.5 mm.Echo tops as measured by high-altitude airborne radar were 7-8 km.The precipitation was driven by a shortwave trough moving from southwest to northeast across the domain.perature (K), relative humidity (%), and winds (m s −1 ) at 00:00 UTC, 31 January, during the middle of the accumulating snowfall.A trough axis is apparent just to the west of the GCPEx domain (green star in Fig. 6).Low-level warmair advection forcing upward motion is coincident with high relative humidity on the leading edge of the trough, over the GCPEx domain (Fig. 6).Temperatures in this layer were around −10 to −15 • C throughout the event, supporting efficient crystal growth, aggregation, and potentially less dense snowfall as this is in the dendritic crystal temperature zone (e.g., Magono and Lee, 1966).Aircraft probe data during a descent over the CARE site between 23:15 and 23:43 UTC showed the median volume diameter (D 0 ) of 3 mm, with particles up to a maximum of 8 mm (aggregates of dendrites) at 2.2 km m.s.l. with a large concentration of smaller sizes < 0.5 mm (dendritic and irregular shapes; Skofronick-Jackson et al., 2015).At the surface, photographs of the precipitation types by the University of Manitoba showed small irregular particles and aggregates (< 3 mm) at 23:30 UTC on 30 January.

D3R radar data
The D3R is a Ku-and Ka-band dual-wavelength polarimetric scanning radar.It was designed for ground validation of rain and falling snow from GPM satellite-borne DPR (dualfrequency precipitation radar).The two frequencies used in the D3R are 13.91 GHz (Ku) and 35.56 GHz (Ka).These two frequencies were used for scattering computations in this research as well.Some parameters of the D3R radar relevant for this paper are shown in Table 2.The range resolution of the radar is adjustable but usually set to 150 m and the nearfield distance is ∼ 300 m; the practical minimum operational range is around 450 m.The minimum detectable signal of the D3R is −10 dBZ at 15 km.This means that when Z h is −10 dBZ at 15 km, the signal-to-noise ratio (SNR) is 0 dB.Therefore, the SNR at any range, r, can be computed as follows: The SNR is a very important indicator for radar data quality control (QC), the other important parameter for QC (in terms of detecting "meteo" vs. "nonmeteo" echoes) being the texture of the standard deviation (SD) of the differential propagation phase (ϕ dp ).We randomly selected 20 out of 85 RHI sweeps from 31 January 2012 and computed the SD of Kuband ϕ dp for each beam over 10 consecutive gates where SNR ≥ 10 dB.According to the histogram of the SD of ϕ dp , 90 % of the values were less than around 8 • .Radar data at a range gate m are identified as "good" data (i.e., meteorological echoes) only if the standard deviation of ϕ dp from the (m − 5)th gate to the (m + 4)th gate is less than 8 • .This criterion sets a good data mask for each beam at Ku-band.
On the other hand, the ϕ dp at Ka-band was determined to be too noisy and hence not used herein.The good data mask for the Ka-band beam is set by the mask determined by the Ku-band criteria, with the additional requirement that the Kaband SNR > 3 dB for the range gate to be considered good.Note that both radars are mounted on a common pedestal so that the Ku and Ka-band beams are perfectly aligned.
There are four scan types that can be performed by the D3R, namely PPI (plan position indicator), RHI (range height indicator), surveillance, and vertical pointing.Figure 7 shows the scan strategies of the D3R on 31 January 2012, which consisted of a fast PPI scan (surveillance scan; 10 • per second) followed by four RHI scans (1 • per second), except from 01:00 to 02:00 UTC.The RHI scans with an azimuth angle of 139.9 • point to the Steam Show site and those at 87.8 • point to the Sky Dive site.There were no RHI scans pointing to the Bob Morton site, and Huronia (52 km) was beyond the operational range (maximum 30 km) of the D3R.During the most intense snowfall the D3R scans did not cover the instrument clusters at the Sky Dive and Steam Show sites.So we were left with the analysis of the D3R radar data at close proximity to the 2DVD or effectively vertical-pointing equivalent using RHI data from 75 to 90 • at the nearest practical range of 600 m.PPI scan data at low elevation angle (3 • ) were also used from range gate at 600 m.The assumption is that there is little evolution of particle microphysics from about 600 m height to the surface and that the synoptic-scale snowfall was uniform in azimuth (confirmed by Skofronick-Jackson et al., 2015).The snowfall was spatially uniform around the CARE site so we selected data at 600 m range to be compared with the 2DVD and Pluvio observations (this range was selected based on the minimum operational range of 450 m; see Table 2) to which 150 m was added based on close examination of data quality.For RHI scans, the Z h at each band was averaged over the beams from 75 to 90 • .The 75 • is obtained from 600•cos(75 • ) ≈ 155 m which is close to the range resolution.For the fast PPI scan, Z h was averaged over all azimuthal beams at 600 m range.
Figure 8 shows the time profile of the averaged Z e at Kuand Ka-bands.There are two problems indicated in this figure.First, theoretically, the Ku-band Z e should be greater than or equal to the Ka-band Z e .The smaller Ku-band Z h indicates that a Z offset exists at both bands.The other problem is that, compared with the Ka-band, there are many dips in the Ku-band Z h .By comparing Fig. 8 with Fig. 7, we found that these dips occur only at RHI scans with azimuth angle larger than 300 • .We examined those RHI scans beam by beam from 90 to 75 • .We further found that when the el-  evation angle is smaller than 78 • , the unreasonably low Z h disappears.Therefore, the RHI scans with azimuth angles larger than 300 • were averaged over the 75 to 78 • elevation angles.To compute the DWR, we need to know the Z offset between the two bands.The measured Z h includes three components (neglecting attenuation): where "error" refers to measurement fluctuations (typically with standard deviation of ∼ 1 dB).The DWR is obtained as the difference between Ku-band Z h and Ka-band Z h , with Z h being in units of dBZ.The measured DWR is as follows: where error (DWR) is now increased, since the Ku-and Kaband measurement fluctuations are uncorrelated (standard deviation of around 1.4 dB).The Z offset is determined by selecting data where the scatterers (snow particles) are sufficiently small in size so that Rayleigh scattering is satisfied at both bands, i.e., DWR true = 0 dB.The criteria are used here to select gates where Ku-band Z h < 0 dBZ along with spatial averaging, which reduces the measurement fluctuations in DWR to estimate Z offset in Eq. (3). Figure 9 shows the averaged Z h for the two bands from 20 RHI scans which satisfy the conditions above.After removing three extreme values (outliers) from Fig. 9, Z offset was estimated as −1.5 dB, which is used in the subsequent data processing.

2DVD data analysis
The 2DVD used in this study was also located at the CARE site.The particle-by-particle mass estimation is based on three methods as follows: 1. Following the procedure in Huang et al. (2015) we use 2DVD single-camera data and apply the weighted Hanesch-matching algorithm (Hanesch, 1999) to rematch snowflakes.A PSD adjustment factor is computed as in Huang et al. (2015) without using the Pluvio gauge as a constraint.Mass is computed from fall speed, D app and environmental conditions using Böhm (1989).The apparent density of the snow (ρ) is defined as 6m/π D 3 app .A mean power-law relation of the form ρ = αD β app is derived for the entire event as in Huang et al. (2015) as well as 1 min averaged N (D app ) is calculated.Note that the scattering model is based on the soft spheroid model with fixed axis ratio = 0.8 and apparent density ρ.The results obtained by this method are denoted HB in the figures and in the rest of the paper.
2. Use the manufacturer's (Joanneum Research, Graz, Austria) matching algorithm and filter-mismatched snowflakes as described in Sect.2.2.The mass is computed from Böhm's equations.The PSD adjustment factor is based on using the Pluvio gauge accumulation as a constraint.Following Liao et al. (2013) as far as the scattering model is concerned, the density is fixed at 0.2 g cc −1 and the volume is computed from mass = density • volume.The effective equal-volume diameter is D eff and the corresponding PSD is denoted N (D eff ), which is different from N (D app ) in ( 1) above.Henceforth, this method is denoted LM.
3. Use Joanneum matching and filtering method as in (2) but compute mass using Heymsfield-Westbrook equations as well as the revised D eff and N (D eff ).This method is denoted HW.Thus, the only difference with ( 2) is in the estimation of mass and the difference in D eff and N (D eff ).The PSD adjustment factor is based on using the Pluvio gauge accumulation as a constraint.The scattering model follows Liao et al. (2013).
The 2DVD measured liquid equivalent snow rate (SR) can be computed directly from mass as follows: where t is the integral time (typically 60 s), N is the number of size bins (typically 101 for the 2DVD), M is the number of snowflakes in the ith size bin, and A j is the measured area of the j th snowflake.Further, V j is the liquid equivalent volume of the j th snowflake, so it is directly related to the mass.Figure 10 compares the liquid equivalent accumulation computed using the three methods above based on 2DVD measurements with the accumulation directly measured by the collocated Pluvio snow gauge.The Pluvio-based accumulation at the end of the event (03:30 Z) was 1.9 mm while the 2DVD-measured accumulations using the three methods are 1.27 mm (HB), 1.45 mm (LM), and 1.24 mm (HW).It is expected that the PSDs of LM and HW should be underestimated because of eliminating mismatched particles which, in principle, could be rematched.Rematching mismatched particles is a research topic on its own and is beyond the scope of this paper.We used a simple way to adjust the PSD for methods 2 and 3 by scaling the PSD by a constant so that the final accumulation matches the Pluvio gauge accumulation.Specifically, the PSD adjustment factors are 1.3 for LM and 1.52 for HW.Note that PSD adjustment of HB (method 1) is not done by forcing 2DVD accumulation to agree with Pluvio, rather the method described in Huang et al. ( 2015) is used giving the adjustment factor of 1.54 for 00:00-00:45 UTC and 1.11 for 00:45-04:00 UTC.From the Pluvio accumulation data in Fig. 10 the SR is nearly constant at 0.7 mm h −1 between relative times of 1.5 and 3 h (or actual time from 01:00 on 30 January to 02:30 UTC).
The radar reflectivities at the two bands are simulated by using the T-matrix method assuming a spheroid shape with an axis ratio of 0.8, consistently with Falconi et al. (2018).The PSD is adjusted for methods 1, 2, and 3 as described above.The orientation angle distribution is assumed to be quasi-random with Gaussian distribution for the zenith angle ] and uniform distribution for the azimuth angle.However, other studies have assumed σ = 10 • (Falconi et al., 2018).The recent observations of snowflake orientation by Garrett et al. (2015) indicate that substantial broadening of the snow orientation distribution can occur due to turbulence. Figure 11 compares the time series of D3Rmeasured Z h with the 2DVD-derived Z e for the entire event (20:00-03:30 UTC at (a) Ku-and (b) Ka-band).The Z e for both bands computed by the three methods generally agree with the D3R measurements to within 3-4 dB.Overall, LM gives the highest Z e and HB gives the lowest, which is especially evident at Ka-band.This is consistent with scattering calculations by Kuo et al. (2016) of single spherical snow aggregates using constant density (0.3 g cc −1 ) giving higher radar cross sections and size-dependent density, i.e., density falls off as inverse size (giving lower cross sections).This feature is consistent with the scattering models referred to herein as LM and HB.
Figure 11.Comparison of the 2DVD-derived Z h with D3R measurements for the entire event, Ku-band (a), and Ka-band (b).This synoptic system started at around 21:00 Z on 30 January and ended at around 03:30 Z on 31 January 2012.Z e by LM is close to HW and slightly higher, whereas the HB method gives the lowest Z e .Z e results computed by all methods generally agree with D3R measured Z h .
From 00:45 to 01:30 UTC on 31 January 2012, the three 2DVD-derived Z e simulations deviate systematically from the D3R results for both bands.The other period is from 23:00 to 23:30 UTC on 30 January 2012, when the Ku-band Z e has significant deviation from the D3R observations but the Ka-band Z e generally agrees with the D3R.Note that this synoptic event started at around 21:00 UTC on 30 January and stopped at 03:30 UTC.We checked the D3R data and found that, before 22:30 UTC, the RHI scans were from 0 to 60 • , so there were no usable data available for comparison with the 2DVD and Pluvio at the CARE site.We note that at 00:30 UTC the King City C-band radar recorded Z h in the range 15-20 dBZ around the CARE site, which is in reasonable agreement with the D3R radar observations (Skofronick-Jackson et al., 2015).
Figure 12a compares the time series of DWR simulated from 2DVD observations with the D3R measurements, whereas Fig. 12b shows the scatterplot In general, HB appears in qualitatively better agreement (better correlated and with significantly less bias) with D3R measurements relative to both LM and HW (significant underestimation relative to D3R).The scatterplot in Fig. 12b is an important result since in the HB method the soft spheroid scattering model is used with density varying approximately inverse with D app (density-D app power law where the larger snow particles have lower density).Hence for a given mass the D app is larger (relative to Ka-band wavelength) and enters the Mie regime, which lowers the radar cross section at Kaband (relative to same mass but constant density radar cross section in LM and HW).Whereas at Ku-band the difference in radar cross sections is less between the two methods (Rayleigh regime).The significant DWR bias in LM and HW relative to DWR observations is somewhat puzzling in that the Liao et al. (2013) scattering model radar cross sections agree with the synthetic complex shaped snow aggregates of the same mass at Ka-band, whereas the HB model underestimates the radar cross section relative to the synthetic complex shaped aggregates.On the other hand, Falconi et al. (2018) demonstrate that the soft spheroid model is adequate at X (close to Ku-band) and Ka-band and by inference adequate for DWR calculations with the caveat that different effective axis ratios may need to be used at Ka-and W-bands.
We also refer to airborne (Ku, Ka) band radar data at 00:30 UTC which showed DWR measurements of 3-6 dB about 1 km height MSL around the CARE site but nearly 0 dB above that all the way to the echo top (Skofronick-Jackson et al., 2015).The latter is not consistent with aircraft spirals over the CARE site about an hour earlier where maximum snow sizes reach ∼ 8 mm.In spite of the difficulty in reconciling the observations from the different sensors, the appropriate scattering model in this particular event appears to favor the soft spheroid model used in HB based on better agreement with DWR observations.The other factor to be considered is the PSD adjustment factor, which is assumed constant and independent of size, which may not be the case, especially for the LM and HW methods as considerable filtering is involved due to mismatch (as discussed in Sect.2.2).Note that a constant PSD adjustment factor will not affect DWR but it will affect Z e .For the HB method Huang et al. (2015) determined the PSD adjustment factor for four events by comparing the 2DVD PSD to that measured by a collocated SVI (snow video imager which was assumed to be the "truth") for each size bin.The PSD adjustment was found to not be size dependent for the HB method.On the other hand, because of the filtering of mismatched particles by the LM and HW methods, the PSD adjustment factor may be size dependent in which case the DWR will also change.More case studies are clearly needed to understand the applicability of the LM and HW methods of simulating DWR.

Snow rate estimation
To obtain radar-SR relationships, we use the 2DVD data and simulations.Since we employ a constant PSD adjustment factor, it will scale both Z e and SR similarly.Figure 13 shows the scatterplot of the 2DVD-derived Z e vs. 2DVD-measured SR along with a power-law fit as Z = aSR b .The fitting method used is based on weighted total least square (WTLS) so the power law can be inverted without any change.The coefficients and exponents of the power-law Z-SR relationship for both bands and three methods are given in Table 3.It is obvious from Fig. 13 that there is considerable scatter at Ku-band for all three methods with the normalized standard deviation (NSD), ranging from 55 % to 70 %.Whereas  of SR is to use Eq. ( 6) when DWR > 1 and SR > 0.2 mm h −1 , else we use the Z Ka -SR power law (note that we do not use the Z Ku -SR power law as the measurement errors of Z Ku seem to be on the high side, Fig. 9).The precise thresholds used herein are ad hoc and may need to be optimized using a much larger data set.Figure 15a shows the radar-derived accumulation using Z Ka -SR vs. the Pluvio accumulation vs. time.The total accumulation from the Pluvio is 2.5 mm and the three radar-based total accumulations, for HB, LM, and HW methods amount to [2.6, 1.8, 2.6 mm].Except for the underestimation in the LM method (−28 %), the other two methods agree with the Pluvio accumulation in this event.
Figure 15b is the same as Fig. 15a, except the combination algorithm mentioned above is used.For this case, ∼ 33 % of data used the Z Ka -SR power law due to threshold constraints given above.The event accumulations for HB, LM, and HW methods amount to [2.4, 1.9, 2.2 mm], which are consistent with the algorithm that uses only the Z Ka -SR power law.However, the criteria of relative bias error in the total accumulation (in events with low accumulations such as this one) are not necessarily an indication that the DWR-based algorithm is not adding value.Rather, the criteria should be snow rate intercomparison, which could not be done due to the low resolution (0.01 mm min −1 ) of the Pluvio 2 400 gauge along with the low event total accumulation of only 2.5 mm.A close qualitative examination of Fig. 15b shows that the HB method more closely "follows" the gauge accumulation relative to HB in time in Fig. 15a.In Fig. 15, the time grid is different for the radar-based data and the gauge data.It is common to linearly interpolate the gauge data to the radar sampling time and if this is done, the rms error for the HB method reduces from 0.1 mm (when using only the Z Ka -SR power law) to 0.045 mm for the DWR algorithm, which constitutes a significant reduction by a factor of 2. The total error in the radar estimate of SR is composed of both parameterization errors as well as measurement errors with measurement errors dominating, since the DWR involves the ratio of two uncorrelated variables.From Sect.8.3 of Bringi and Chandrasekar (2001), the total error of SR in Eq. ( 6) is around 50 % (ratio of standard deviation to the mean).The assumptions are (a) the standard deviation of the measurement of Z e is 0.8 dB, (b) the standard deviation of the DWR (in dB) measurement is 1.13 dB, and (c) the parameterization error is 30 % from Table 4.However, considering the The radar SR is computed by combining SR(Z Ku ,DWR) and Z Ka -SR as described in the text.The accumulated SR derived from the radar using HB method is 2.38 mm, using LM it is 1.94 mm and using HW it is 2.24 mm.Z e fluctuations in Fig. 9, the measurement standard deviation probably exceeds 0.8 dB, especially at Ku-band.Thus, sufficient smoothing of DWR is needed to minimize the measurement error as much as possible while maintaining sufficient spatial resolution.
Note that the error model used here is additive with the parameterization, and measurement errors modeled as zero mean and uncorrelated with the corresponding error variances estimated either from data or via simulations (as described in Sect.7 of Bringi and Chandrasekar, 2001).This is a simplified error model since it assumes that radar Z and snow gage measurements are unbiased based on accurate calibration.A more elaborate approach of quantifying uncertainty in precipitation rates is described by Kirstetter et al. (2015).

Summary and conclusions
The main objective of this paper is to develop a technique for snow estimation using scanning dual-wavelength radar operating at Ku-and Ka-bands (D3R radar operated by NASA).We use the 2-D video disdrometer and collocated Pluvio gauge to derive an algorithm to retrieve snow rate from reflectivity measurements at the two frequencies compared to the conventional single-frequency Z e -SR power laws.The important microphysical information needed is provided by the 2DVD to estimate the mass of each particle knowing the fall speed, apparent volume, area ratio, and environmental factors from which an average density-size relation is derived (e.g., Huang et al., 2015;von Lerber et al., 2017;Böhm, 1989;Heymsfield and Westbrook, 2010).
We describe in detail the data processing of 2DVD camera images (in two orthogonal planes) and the role of particle mismatches that give erroneous fall speeds.We use the Huang et al. (2015) method of rematching using singlecamera data but also use the manufacturer's matching code with substantial filtering of the mismatched particles since the apparent volume and diameter (D app ) are more accurate.To account for the filtering of the mismatched particles, the particle size distribution (in methods 2 and 3 in Sect.3.3) is adjusted by a constant factor using the total accumulation from the Pluvio as a constraint.
Two scattering models are used to compute the Z Ku and Z Ka , termed the soft spheroid model (Huang et al., 2015;HB method) and the Liao-Meneghini (LM) model, which uses the concept of effective density.In these two methods the particle mass is based on Böhm (1989).The method of Heymsfield and Westbrook (2010) is also used to estimate mass which is similar to Böhm (1989) but is expected to be more accurate (Westbrook and Sephton, 2017); along with the LM model for scattering, this method is termed HW.
The case study chosen is a large-scale synoptic snow event that occurred over the instrumented site of CARE during GCPEx.The Z Ku and Z Ka were simulated based on 2DVD data and the three methods, i.e., HB, LM, and HW yielded similar values within ±3 dB.When compared with D3R radar measurements extracted as a time series over the instrumented site, the LM and HW methods were closer to the radar measurements with the HB method being lower by ≈ 3 dB.Some systematic deviations of simulated reflectivities by the three methods from the radar measurements were explained by a possible size dependence of the PSD adjustment factor.
The direct comparison of DWR (ratio of Z Ku to Z Ka ) from simulations with DWR measured by radar showed that the HB method gave the lowest bias with the data points more or less evenly distributed along the 1 : 1 line.The simulation of DWR by LM and HW methods underestimated the radar measurements of DWR quite substantially, even though the correlation appeared to be reasonable.The reason for this discrepancy is difficult to explain since a constant PSD adjustment factor (different for method 1 relative to methods 2 and 3 in Sect.3.3) would not affect the DWR.From the scattering model viewpoint, the LM method takes into account the complex shapes of snow aggregates via an effective density approach, whereas the HB method uses soft spheroid model with density varying approximately inversely with size.We did not attempt to classify the particle types in this study.
The retrieval of SR was formulated as SR = c•Z d Ku •DWR e , where [c, d, e] were obtained via nonlinear least squares for the three methods.The total accumulation from the three methods using radar-measured Z Ku and DWR were compared with the total accumulation from the Pluvio (2.5 mm) to demonstrate closure.The closest to Pluvio was the HB method (2.4 mm), next was the HW method (2.24 mm) and then there was LM (1.94 mm).At such low total accumulations, the three methods show good agreement with each other as well as with the Pluvio gauge.The poor resolution of the gauge combined with the relatively low total accumulation in this event precluded direct comparison of snow rates.The combined estimate of parameterization and measurement errors for snow rate estimation was around 50 %.From variance decomposition, the measurement error variance as a fraction of the total error variance was 58 %, and the parameterization error variance fraction was 42 %.Further, the DWR was responsible for 90 % of the measurement error variance, which is not surprising since it is the ratio of two uncorrelated reflectivities.Thus, the DWR radar data have to be smoothed spatially (in range and azimuth) to reduce this error, which will degrade the spatial resolution but is not expected to pose a problem in large-scale synoptic snow events.
The snow rate estimation algorithms developed here are expected to be applicable to similar synoptic-forced snowfall under similar environmental conditions (e.g., temperature and relative humidity) but not, for example, to lake effect snowfall as the microphysics are quite different.However, analyses of more events are needed before any firm conclusions can be drawn as to applicability to other regions or environmental conditions.
Data availability.The data used in this study can be made available upon request to the corresponding author.

Figure 1 .
Figure 1.A snowflake observed by a 2DVD from two views.The thick black line is the contour of the snowflake and the thin black lines show the holes inside the snowflake.The effective area, A e , equals the area enclosed by the thick black curve minus the area enclosed by thin lines.The blue line represents the minimum circumscribed ellipse, the enclosed area of which is denoted by A.

Figure 3 .
Figure 3. Fall speed vs. D app for the synoptic case on 31 January 2012 at the CARE site.The green circles represent the results of the manufacturer's matching algorithm, which is known to allow mismatched particles with unrealistic fall speeds.The first filtering step is the selection of matched particles which satisfy Hanesch scan line criteria (magenta).The second filter step is shown as blue crosses, which are based on particles with density (from mass computed by Böhm's or Heymsfield-Westbrook method) lower than 1 g cm −3 .

Figure 4 .
Figure 4.A map of the GCPEx field campaign.The five test sites are CARE, Sky Dive, Steam Show, Bob Morton, and Huronia.The ground observation instruments, namely 2DVD, D3R, and Pluvio, used in this research, were located at CARE.

Figure 6 .
Figure 6.The 00:00 UTC, 31 January 2012, 850 hPa geopotential heights (m, black solid contours), temperature (K, red: above freezing, blue: below freezing), relative humidity (%, green shaded contours), and wind (m s −1 , wind barbs).The red dot in the center right portion of the figure denotes the general location of the GCPEx field instruments.

Figure 8 .
Figure 8.The time series of averaged raw Z h at the CARE site.There are two problems indicated in this figure: (i) the Ku-band Z h is smaller than the Ka-band Z h on average.(ii) Compared with the Ka-band, there are many too small values of the Ku-band Z h .

Figure 9 .
Figure 9.The averaged raw Z h for Ku-and Ka-bands.The Z h was randomly selected from 20 of 85 RHI scans with Ku-band Z h < 0 dBZ, range < 1 km, and Ka-band SNR > 3 dB.

Figure 10 .
Figure10.Comparison of liquid equivalent accumulations computed using HB, LM, and HW methods based on 2DVD measurements and that were directly measured by the collocated Pluvio snow gauge.We used the total accumulation to estimate the PSD adjustment factor for the LM and HW methods.

Figure 12 .
Figure 12.Comparison of the 2DVD-derived DWR using HB, LM, and HW methods with the D3R-measured DWR.Panel (a) shows the time profile of the D3R, and (b) shows the scatterplot.

Figure 14 .
Figure 14.Estimated SR using Z e and DWR of the 2DVD and Eq.(10) vs. 2DVD SR scatterplot for the HB method (a), LM method (b), and HW method (c).

Figure 15 .
Figure 15.Comparison of the radar-derived accumulated SR using HB, LM, and HW methods with Pluvio gauge measurement.(a) The radar SR is computed by Z Ka -SR relationships.The Pluvioaccumulated SR on 03:18 UTC is 2.48 mm.The radar-accumulated SRs for HB, LM, and HW are 2.64, 1.81, and 2.66 mm.(b)The radar SR is computed by combining SR(Z Ku ,DWR) and Z Ka -SR as described in the text.The accumulated SR derived from the radar using HB method is 2.38 mm, using LM it is 1.94 mm and using HW it is 2.24 mm.

Table 4 .
Coefficients and exponents of the SR(Z Ku ,DWR) relation (see Eq. 10) for three methods.

Table A1 .
radar technique development for snow rate estimation Appendix A List of acronyms and symbols.Equivalent reflectivity factor.In this paper, we refer to the reflectivity factor based on scattering computation.