Research article 01 Mar 2019
Research article  01 Mar 2019
Dualwavelength radar technique development for snow rate estimation: a case study from GCPEx
 ^{1}Center for Atmospheric REmote sensing (CARE), Kyungpook National University, Daegu 41566, Republic of Korea
 ^{2}Department of Electrical and Computer Engineering, Colorado State University, Fort Collins, Colorado 80523, USA
 ^{3}Research Applications Lab, National Center for Atmospheric Research, Boulder, Colorado 80307, USA
 ^{4}Institute for Atmospheric and Earth System Research, University of Helsinki, 00560 Helsinki, Finland
 ^{5}Finnish Meteorological Institute, 00101 Helsinki, Finland
 ^{1}Center for Atmospheric REmote sensing (CARE), Kyungpook National University, Daegu 41566, Republic of Korea
 ^{2}Department of Electrical and Computer Engineering, Colorado State University, Fort Collins, Colorado 80523, USA
 ^{3}Research Applications Lab, National Center for Atmospheric Research, Boulder, Colorado 80307, USA
 ^{4}Institute for Atmospheric and Earth System Research, University of Helsinki, 00560 Helsinki, Finland
 ^{5}Finnish Meteorological Institute, 00101 Helsinki, Finland
Correspondence: GwoJong Huang (gwojong.huang@colostate.edu)
Hide author detailsCorrespondence: GwoJong Huang (gwojong.huang@colostate.edu)
quantitative precipitation estimation (QPE) of snowfall has generally been expressed in powerlaw form between equivalent radar reflectivity factor (Z_{e}) and liquid equivalent snow rate (SR). It is known that there is large variability in the prefactor 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 dualwavelength radar reflectivity ratio (DWR) technique can improve SR accuracy by estimating one of the PSD parameters (characteristic diameter), thus reducing the variability due to the prefactor. The two frequencies commonly used in dualwavelength techniques are Ku and Kabands. The basic idea of DWR is that the snow particle sizetowavelength ratio is falls in the Rayleigh region at Kuband but in the Mie region at Kaband.
We propose a method for snow rate estimation by using NASA D3R radar DWR and Kaband reflectivity observations collected during a longduration synoptic snow event on 30–31 January 2012 during the GCPEx (GPM Coldseason 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 it directly with D3R radar measurements. We also use the 2DVDbased mass to compute the 2DVDbased SR. Using three different mass estimation methods, we arrive at three respective sets of Z–SR and SR(Z_{h}, DWR) relationships. We then use these relationships with D3R measurements to compute radarbased SR. Finally, we validate our method by comparing the D3R radarretrieved SR with accumulated SR directly measured by a wellshielded Pluvio gauge for the entire synoptic event.
A detailed understanding of the geometric, microphysical, and scattering properties of ice hydrometeors is a vital prerequisite for the development of radarbased quantitative precipitation estimation (QPE) algorithms. Recent advances in surface and airborne optical imaging instruments and the wide proliferation of dualpolarization and multiwavelength radar systems (ground based, airborne or satellite) have allowed for observations of the complexity inherent in winter precipitation via dedicated field programs (e.g., SkofronickJackson 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 coexistence 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) 2D video disdrometer (2DVD; Schönhuber et al., 2008), (ii) precipitation imaging package (PIP; von Lerber et al., 2017), (iii) MultiAngle Snowflake Camera (MASC; Garrett et al., 2012). When these instruments are used in conjunction with a wellshielded GEONOR or PLUVIO gauge, it is shown that a physically consistent representation of the geometric, microphysical, and scattering properties needed for radarbased 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.
Radarbased 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 Wband) 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 Wbands) (e.g., Heymsfield et al., 2005, 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 dualwavelength reflectivity ratio (DWR, the ratio of reflectivity from two different bands) radarbased 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 verticalpointing ground radars or nadirpointing airborne radars (Liao et al., 2005, 2008, 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 sizetowavelength ratio falls in the Rayleigh region at a lowfrequency band (e.g., Kuband) but in the Mie region at a highfrequency band (e.g., Kaband) (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 liquidequivalent size or mass (Liao et al., 2016). In this sense the DWR is similar to differential reflectivity (Z_{dr}) in dualpolarization 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 dualwavelength 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 DWRbased 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.
2.1 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 powerlaw 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 methodology has been described in detail by Szyrmer and Zawadzki (2010), Huang et al. (2015), and von 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}_{\mathrm{r}}^{\mathrm{0.25}}$ (${A}_{\mathrm{r}}={A}_{\mathrm{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}_{\mathrm{r}}^{\mathrm{0.5}}$ 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 shapecorrecting 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.
2.2 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 twodimensional 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 linescan cameras (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 groundbased instruments that can automatically and continuously capture the horizontal projected views (i.e., in the plane orthogonal to the flow) of precipitation particles (however, 3Dreconstruction based on multiple views can give this information; Kleinkort et al., 2017). Compared with other opticalbased 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}, respectively, the distance between two optic planes is D_{d}, and camera A observed a particle at t_{0}, we have ${t}_{\mathrm{1}}={t}_{\mathrm{0}}+{D}_{\mathrm{d}}/{v}_{\mathrm{max}}$ and ${t}_{\mathrm{2}}={t}_{\mathrm{0}}+{D}_{\mathrm{d}}/{v}_{\mathrm{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 redid 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 singlecamera 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 diameter 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 ($=\mathit{\pi}{D}_{\mathrm{app}}^{\mathrm{3}}/\mathrm{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 twostep 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.
2.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 Tmatrix 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 axissymmetric 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. (2006, 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 quasirandom 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 Wband (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 Tmatrix to compute Z_{e} (at X, Ka, and Wbands) 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 Xband, 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 complexshaped aggregates to the soft spheroid model at Wband 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. complexshape 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 Kabands, 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 equalmass spheroid, Liao et al. (2013) used sixbranch 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 (Kaband), 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 complexshaped 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 Kabands (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 (samemass) complexshaped aggregate snowflake (Liao et al., 2016).
3.1 Test site instrumentation and the synoptic event
The GPM Coldseason 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 multifrequency active and passive microwave sensors to detect and estimate falling snow…” (SkofronickJackson et al., 2015). The field experiment sites were located north of Toronto, Canada between Lake Huron and Lake Ontario. The GCPEx had five test sites, namely CARE (Centre for Atmospheric Research Experiments), Sky Dive, Steam Show, Bob Morton, and Huronia. The locations of five sites are shown in Fig. 4. The CARE site was the main test site for the experiment, located at 44^{∘}13^{′}58.44^{′′} N, 79^{∘}46^{′}53.28^{′′} W and equipped with an extensive suite of ground instruments. The 2DVD (SN37) and OTT Pluvio^{2} 400 used for observations and analyses in this paper were installed inside a DFIR (double fence intercomparison reference) wind shield. The dualfrequency dualpolarized doppler radar (D3R) was also located at the CARE site (Vega et al., 2014) near the 2DVD. The instruments used in this paper are depicted in Fig. 5. Because the radar and the instrumented site were nearly collocated, we can effectively view the setup as similar to a verticalpointing radar as described in more detail in Sect. 3.2.
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 SkofronickJackson 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 highaltitude airborne radar were 7–8 km. The precipitation was driven by a shortwave trough moving from southwest to northeast across the domain. Figure 6 displays the 850 hPa geopotential heights (m), temperature (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). Lowlevel 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; SkofronickJackson 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.
3.2 D3R radar data
The D3R is a Ku and Kaband dualwavelength polarimetric scanning radar. It was designed for ground validation of rain and falling snow from GPM satelliteborne 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 signaltonoise 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 Kuband. On the other hand, the φ_{dp} at Kaband was determined to be too noisy and hence not used herein. The good data mask for the Kaband beam is set by the mask determined by the Kuband 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 Kaband 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 verticalpointing 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 synopticscale snowfall was uniform in azimuth (confirmed by SkofronickJackson 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 $\mathrm{600}\cdot \mathrm{cos}\left(\mathrm{75}{}^{\circ}\right)\approx \mathrm{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 Ku and Kabands. There are two problems indicated in this figure. First, theoretically, the Kuband Z_{e} should be greater than or equal to the Kaband Z_{e}. The smaller Kuband Z_{h} indicates that a Z offset exists at both bands. The other problem is that, compared with the Kaband, there are many dips in the Kuband 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 elevation 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 Kuband Z_{h} and Kaband 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 Kuband 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.
3.3 2DVD data analysis
The 2DVD used in this study was also located at the CARE site. The particlebyparticle mass estimation is based on three methods as follows:

Following the procedure in Huang et al. (2015) we use 2DVD singlecamera data and apply the weighted Haneschmatching 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 $\mathrm{6}m/\mathit{\pi}{D}_{\mathrm{app}}^{\mathrm{3}}$. A mean powerlaw relation of the form $\mathit{\rho}=\mathit{\alpha}{D}_{\mathrm{app}}^{\mathit{\beta}}$ 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.

Use the manufacturer's (Joanneum Research, Graz, Austria) matching algorithm and filtermismatched 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 $\text{mass}=\text{density}\cdot \text{volume}$. The effective equalvolume 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.

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 jth snowflake. Further, V_{j} is the liquid equivalent volume of the jth 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 Pluviobased accumulation at the end of the event (03:30 Z) was 1.9 mm while the 2DVDmeasured 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 Tmatrix 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 quasirandom with Gaussian distribution for the zenith angle $[\text{mean}=\mathrm{0}{}^{\circ},\phantom{\rule{0.25em}{0ex}}\mathit{\sigma}=\mathrm{45}{}^{\circ}]$ and uniform distribution for the azimuth angle. However, other studies have assumed $\mathit{\sigma}=\mathrm{10}{}^{\circ}$ (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 2DVDderived Z_{e} for the entire event (20:00–03:30 UTC at (a) Ku and (b) Kaband). 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 Kaband. 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 sizedependent 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.
From 00:45 to 01:30 UTC on 31 January 2012, the three 2DVDderived 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 Kuband Z_{e} has significant deviation from the D3R observations but the Kaband 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 Cband radar recorded Z_{h} in the range 15–20 dBZ around the CARE site, which is in reasonable agreement with the D3R radar observations (SkofronickJackson 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} (densityD_{app} power law where the larger snow particles have lower density). Hence for a given mass the D_{app} is larger (relative to Kaband 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 Kuband 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 Kaband, 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 Kuband) and Kaband and by inference adequate for DWR calculations with the caveat that different effective axis ratios may need to be used at Ka and Wbands.
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 (SkofronickJackson 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.
3.4 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 2DVDderived Z_{e} vs. 2DVDmeasured SR along with a powerlaw 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 powerlaw 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 Kuband for all three methods with the normalized standard deviation (NSD), ranging from 55 % to 70 %. Whereas at Kaband the scatter is significantly lower with NSD from 40 % to 45 %. The errors in Table 3 are generally termed parameterization errors.
By using dualwavelength radar, we can estimate SR using Z_{e} at two bands as follows:
where ${a}^{\prime}=(\mathrm{1}/a{)}^{{b}^{\prime}}$ and ${b}^{\prime}=\mathrm{1}/b$. To reduce error, we may take the geometric mean of these two estimators as follows:
where $c=({a}_{\mathrm{1}}^{\prime}{a}_{\mathrm{2}}^{\prime}{)}^{\mathrm{1}/\mathrm{2}}$, $d=({b}_{\mathrm{1}}^{\prime}+{b}_{\mathrm{2}}^{\prime})/\mathrm{2}$, and $e={b}_{\mathrm{2}}^{\prime}/\mathrm{2}$. Note that the DWR in Eq. (6) is on a linear scale, i.e., expressed as a ratio of reflectivity in units of mm^{6} m^{−3}. Using Table 3 to set the initial guess of (c, d, e), nonlinear least squares fitting was used to determine the optimized (c, d, e) with the cost function being the squared difference between the 2DVDbased measurements of SR and $c{Z}_{\mathrm{Ku}}^{d}{\text{DWR}}^{e}$, where Z_{Ku} and DWR are from 2DVD simulations. Figure 14 shows the SR computed from the 2DVD simulations of Kuband Z_{e} and the DWR using Eq. (6) vs. the 2DVDmeasured SR. The (c, d, e) values for the three methods are given in Table 4. As can also be seen from Fig. 14 and Table 4, the SR(Z_{Ku},DWR) using the LM method results in the lowest NSD of 28.49 %, but the other two methods have similar values of NSD (≈30 %) and, as such, these differences are not statistically significant. Although SR(Z_{Ku},DWR) has a smaller parameterization error than Z_{e}–SR, the SR(Z_{Ku},DWR) estimation is biased high when SR<0.2 mm h^{−1} (see Fig. 14). When SR is small, the size of snowflakes is usually also small and falls in the Rayleigh region at both frequencies, resulting in DWR very close to 1 (when expressed as a ratio). This implies that there is no information content in the DWR so including it just adds to the measurement error. Hence, for small SR or when DWR ≈1, we use the Z_{e}–SR power law.
So far the singlefrequency SR retrieval algorithms were based on 2DVDbased simulations with a PSD adjustment factor using the total accumulation from Pluvio as a constraint. The algorithm we propose for radarbased estimation 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 radarderived accumulation using Z_{Ka}–SR vs. the Pluvio accumulation vs. time. The total accumulation from the Pluvio is 2.5 mm and the three radarbased 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 DWRbased 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 radarbased 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 Z_{e} fluctuations in Fig. 9, the measurement standard deviation probably exceeds 0.8 dB, especially at Kuband. 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).
The main objective of this paper is to develop a technique for snow estimation using scanning dualwavelength radar operating at Ku and Kabands (D3R radar operated by NASA). We use the 2D 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 singlefrequency 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 densitysize 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 largescale 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 $\text{SR}=c\cdot {Z}_{\mathrm{Ku}}^{d}\cdot {\text{DWR}}^{e}$, where [$c,d,e$] were obtained via nonlinear least squares for the three methods. The total accumulation from the three methods using radarmeasured 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 largescale synoptic snow events.
The snow rate estimation algorithms developed here are expected to be applicable to similar synopticforced 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.
The data used in this study can be made available upon request to the corresponding author.
GJH and VNB developed the main idea of this article. GJH analyzed 2DVD data, processed D3R data, adapted the scattering models for the case considered, simulated radar reflectivity factor, derived Z_{e}–SR and SR(Z,DWR) relationships, and wrote most of the article draft. AJN described the meteorological aspects of the studied case and wrote Sect. 3.1. VNB, DM, and BMN validated the scattering model used in the article. VNB and GL validated the D3R data. VNB and BMN supervised the NASA PMM Science research. GL supervised Korea Brain Pool Program and KMA grant. BMN and VNB reviewed and edited the article.
The authors declare that they have no conflict of interest.
GwoJong Huang acknowledges support from the Brain Pool Program through the National
Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT
(grant number 171S531874). All authors except Gyuwon Lee acknowledge support
from NASA PMM Science grant NNX16AE43G. Authors GwoJong Huang and Gyuwon Lee were funded by
the Korea Meteorological Administration Research and Development Program
under Grant KMI201806810.
Edited by: S. Joseph Munchak
Reviewed by: three anonymous referees
Abraham, F. F.: Functional dependence of drag coefficient of a sphere on Reynolds number, Phys. Fluids, 13, 2194–2195, https://doi.org/10.1063/1.1693218, 1970.
Barthazy, E., Göke, S., Schefold, R., and Högl, D.: An optical array instrument for shape and fall velocity measurements of hydrometeors, J. Atmos. Ocean. Tech., 21, 1400–1416, https://doi.org/10.1175/JTECHD160221.1, 2004.
Bernauer, F., Hürkamp, K., Rühm, W., and Tschiersch, J.: On the consistency of 2D video disdrometers in measuring microphysical parameters of solid precipitation, Atmos. Meas. Tech., 8, 3251–3261, https://doi.org/10.5194/amt832512015, 2015.
Böhm, H. P.: A general equation for the terminal fall speed of solid hydrometeors, J. Atmos. Sci., 46, 2419–2427, https://doi.org/10.1175/15200469(1989)046<2419:AGEFTT>2.0.CO;2, 1989.
Botta, G., Aydin, K., and Verlinde, J.: Modeling of microwave scattering from cloud ice crystal aggregates and melting aggregates: A new approach, IEEE Geosci. Remote Sens. Lett., 7, 572–576, https://doi.org/10.1109/LGRS.2010.2041633, 2010.
Brandes, E. A., Ikeda, K., Zhang, G., Schönhuber, M., and Rasmussen, R. M.: A statistical and physical description of hydrometeor distributions in Colorado snowstorms using a video disdrometer, J. Appl. Meteor. Climatol., 46, 634–650, https://doi.org/10.1175/JAM2489.1, 2007.
Bringi, V. N. and Chandrasekar, V.: Polarimetric Doppler Weather Radar: Principles and Applications, Cambridge University Press, Cambridge, United Kingdom, 545–554, 2001.
Bukovčić, P., Ryzhkov, A., Zrnić, D., and Zhang, G.: Polarimetric Radar Relations for Quantification of Snow Based on Disdrometer Data, J. Appl. Meteor. Climatol., 57, 103–120, https://doi.org/10.1175/JAMCD170090.1, 2018.
Chobanyan, E., Šekeljić, N. J., Manić, A. B., Ilić, M. M., Bringi, V. N., and Notaroš, B. M.: Efficient and Accurate Computational Electromagnetics Approach to Precipitation Particle Scattering Analysis Based on HigherOrder Method of Moments Integral Equation Modeling, J. Atmos. Ocean. Tech., 32, 1745–1758, https://doi.org/10.1175/JTECHD150037.1, 2015.
Davies, C. N.: Definitive Equations for the Fluid Resistance of Spheres, Proc. Phys. Soc. London, A57, 259–270, https://doi.org/10.1088/09595309/57/4/301, 1945.
Delanoë, J., Protat, A., Testud, J., Bouniol, D., Heymsfield, A. J., Bansemer, A., Brown, P. R. A., and Forbes, R. M.: Statistical properties of the normalized ice particle size distribution, J. Geophys. Res., 110, D10201, https://doi.org/10.1029/2004JD005405, 2005.
Draine, B. T. and Flatau, P. J.: Discrete dipole approximation for scattering calculations, J. Opt. Soc. Am. A, 11, 1491–1499, https://doi.org/10.1364/JOSAA.11.001491, 1994.
Falconi, M. T., von Lerber, A., Ori, D., Marzano, F. S., and Moisseev, D.: Snowfall retrieval at X, Ka and W bands: consistency of backscattering and microphysical properties using BAECC groundbased measurements, Atmos. Meas. Tech., 11, 3059–3079, https://doi.org/10.5194/amt1130592018, 2018.
Fujiyoshi, Y., Endoh, T., Yamada, T., Tsuboki, K., Tachibana, Y., and Wakahama, G.: Determination of a ZR relationship or snowfall using a radar and high sensitivity snow gauges, J. Appl. Meteor., 29, 147–152, https://doi.org/10.1175/15200450(1990)029<0147:DOARFS>2.0.CO;2, 1990.
Garrett, T. J., Fallgatter, C., Shkurko, K., and Howlett, D.: Fall speed measurement and highresolution multiangle photography of hydrometeors in free fall, Atmos. Meas. Tech., 5, 2625–2633, https://doi.org/10.5194/amt526252012, 2012.
Garrett, T. J., Yuter, S. E., Fallgatter, C., Shkurko, K., Rhodes, S. R., and Endries, J. L.: Orientations and Aspect ratios of falling snow, Geophys. Res. Lett., 42, 4617–4622, https://doi.org/10.1002/2015GL064040, 2015.
Hanesch, M.: Fall Velocity and Shape of Snowflake, PhD dissertation, Swiss Federal Institute of Technology, Zürich, Switzerland, available at: https://doi.org/10.3929/ethza003837623 (last access: 19 February 2019), 1999.
Heymsfield, A. J. and Westbrook, C. D.: Advances in the Estimation of Ice Particle Fall Speeds Using Laboratory and Field Measurements, J. Atmos. Sci., 67, 2469–2482, https://doi.org/10.1175/2010JAS3379.1, 2010.
Heymsfield, A. J., Wang, Z., and Matrosov, S.: Improved Radar Ice Water Content Retrieval Algorithms Using Coincident Microphysical and Radar Measurements, J. Appl. Meteor., 44, 1391–1412, https://doi.org/10.1175/JAM2282.1, 2005.
Heymsfield, A. J., Krämer, M., Wood, N. B., Gettelman, A., Field, P. R., and Liu, G.: Dependence of the Ice Water Content and Snowfall Rate on Temperature, Globally: Comparison of in Situ Observations, Satellite Active Remote Sensing Retrievals, and Global Climate Model Simulations, J. Appl. Meteor. Climatol., 56, 189–215, https://doi.org/10.1175/JAMCD160230.1, 2016.
Hogan, R. J., Mittermaier, M. P., and Illingworth, A. J.: The retrieval of ice water content from radar reflectivity factor and temperature and its use in evaluating a mesoscale model, J. Appl. Meteor. Climatol., 45, 301–317, https://doi.org/10.1175/JAM2340.1, 2006.
Huang, G., Bringi, V. N., Cifelli, R., Hudak, D., and Petersen, W. A.: A Methodology to Derive Radar Reflectivity–Liquid Equivalent Snow Rate Relations Using CBand Radar and a 2D Video Disdrometer, J. Atmos. Ocean. Tech., 27, 637–651, https://doi.org/10.1175/2009JTECHA1284.1, 2010.
Huang, G., Bringi, V. N., Moisseev, D., Petersen, W. A., Bliven, L., and Hudak, D.: Use of 2DVideo Disdrometer to Derive Mean DensitySize and Z_{e}SR Relations: Four Snow Cases from the Light Precipitation Validation Experiment, Atmos. Res., 153, 34–48, https://doi.org/10.1016/j.atmosres.2014.07.013, 2015.
Kirstetter, P. E., Gourley, J. J., Hong, Y., Zhang, J., Moazamigoodarzi, S., Langston, C., and Arthur, A.: Probabilistic Precipitation Rate Estimates with Groundbased Radar Networks, Water Resour. Res., 51, 1422–1442, https://doi.org/10.1002/2014WR015672, 2015.
Kleinkort, C., Huang, G., Bringi, V. N., and Notaroš, B. M.: Visual Hull Method for Realistic 3D Particle Shape Reconstruction Based on HighResolution Photographs of Snowflakes in Free Fall from Multiple Views, J. Atmos. Ocean. Tech., 34, 679–702, https://doi.org/10.1175/JTECHD160099.1, 2017.
Kneifel, S., von Lerber, A., Tiira, J., Moisseev, D., Kollias, P., and Leinonen, J.: Observed relations between snowfall microphysics and triplefrequency radar measurements, J. Geophys. Res.Atmos., 120, 6034–6055, https://doi.org/10.1002/2015JD023156, 2015
Kuo, K., Olson, W. S., Johnson, B. T., Grecu, M., Tian, L., Clune, T. L., van Aartsen, B. H., Heymsfield, A. J., Liao, L., and Meneghini, R.: The Microwave Radiative Properties of Falling Snow Derived from Nonspherical Ice Particle Models. Part I: An Extensive Database of Simulated Pristine Crystals and Aggregate Particles, and Their Scattering Properties, J. Appl. Meteor. Climatol., 55, 691–708, https://doi.org/10.1175/JAMCD150130.1, 2016.
Leinonen, J., Kneifel, S., Moisseev, D., Tyynelä, J., Tanelli, S., and Nousiainen, T.: Evidence of nonspheroidal behavior in millimeterwavelength radar observations of snowfall, J. Geophys. Res., 117, D18205, https://doi.org/10.1029/2012JD017680, 2012.
Liao, L., Meneghini, R., Iguchi, T., and Detwiler, A.: Use of DualWavelength Radar for Snow Parameter Estimates, J. Atmos. Ocean. Tech., 22, 1494–1506, https://doi.org/10.1175/JTECH1808.1, 2005.
Liao, L., Meneghini, R., Tian, L., and Heymsfield, G. M.: Retrieval of snow and rain from combined X and Wband airborne radar measurements, IEEE T. Geosci. Remote, 46, 1514–1524, https://doi.org/10.1109/TGRS.2008.916079, 2008.
Liao, L., Meneghini, R., Nowell, H. K., and Liu, G.: Scattering computations of snow aggregates from simple geometrical particle models, IEEE J. Sel. Top. Earth Obs. Remote Sens., 6, 1409–1417, https://doi.org/10.1109/JSTARS.2013.2255262, 2013.
Liao, L., Meneghini, R., Tokay, A., and Bliven, L. F.: Retrieval of snow properties for Ku and Kaband dualfrequency radar, J. Appl. Meteor. Climatol., 55, 1845–1858, https://doi.org/10.1175/JAMCD150355.1, 2016.
Magono, C. and Lee, C. W.: Meteorological Classification of Natural Snow Crystals, J. Fac. Sci., Hokkaidô University, Japan, Ser., 7, 321–335, 1966.
Matrosov, S. Y.: A DualWavelength Radar Method to Measure Snowfall Rate, J. Appl. Meteor., 37, 1510–1521, https://doi.org/10.1175/15200450(1998)037<1510:ADWRMT>2.0.CO;2, 1998.
Matrosov, S. Y., Heymsfield, A. J., and Wang, Z.: Dualfrequency radar ratio of nonspherical atmospheric hydrometeors, Geophys. Res. Lett., 32, L13816, https://doi.org/10.1029/2005GL023210, 2005.
Meneghini, R. and Liao, L.: On the Equivalence of DualWavelength and DualPolarization Equations for Estimation of the Raindrop Size Distribution, J. Atmos. Ocean. Tech., 24, 806–820, https://doi.org/10.1175/JTECH2005.1, 2007.
Mishchenko, M. I., Travis, L. D., and Lacis, A. A.: Scattering, Absorption, and Emission of Light by Small Particles, Cambridge University Press, Cambridge, United Kingdom, 2002.
Newman, A. J., Kucera, P. A., and Bliven, L. F.: Presenting the Snowflake Video Imager (SVI), J. Atmos. Ocean. Tech., 26, 167–179, https://doi.org/10.1175/2008JTECHA1148.1, 2009.
Petäjä, T., O'Connor, E., Moisseev, D., Sinclair, V., Manninen, A., Väänänen, R., von Lerber, A., Thornton, J., Nicoll, K., Petersen, W., Chandrasekar, V., Smith, J., Winkler, P., Krüger, O., Hakola, H., Timonen, H., Brus, D., Laurila, T., Asmi, E., Riekkola, M., Mona, L., Massoli, P., Engelmann, R., Komppula, M., Wang, J., Kuang, C., Bäck, J., Virtanen, A., Levula, J., Ritsche, M., and Hickmon, N.: BAECC: A Field Campaign to Elucidate the Impact of Biogenic Aerosols on Clouds and Climate, B. Am. Meteorol. Soc., 97, 1909–1928, https://doi.org/10.1175/BAMSD1400199.1, 2016.
Petty, G. W. and Huang, W.: Microwave backscattering and extinction by soft ice spheres and complex snow aggregates, J. Atmos. Sci., 67, 769–787, https://doi.org/10.1175/2009JAS3146.1, 2010.
Rasmussen, R., Dixon, M., Vasiloff, S., Hage, F., Knight, S., Vivekanandan, J., and Xu, M.: Snow Nowcasting Using a RealTime Correlation of Radar Reflectivity with Snow Gauge Accumulation, J. Atmos. Sci., 42, 20–36, https://doi.org/10.1175/15200450(2003)042<0020:SNUART>2.0.CO;2, 2003.
Ryzhkov, A. V., Zrnić, D. S., and Gordon, B. A.: Polarimetric Method for Ice Water Content Determination, J. Appl. Meteor., 37, 125–134, https://doi.org/10.1175/15200450(1998)037<0125:PMFIWC>2.0.CO;2, 1998.
Schönhuber, M., Urban, H. E., Randeu, W. L., and Poiares Baptista, J. P. V.: Empirical Relationships between Shape, Water Content and Fall Velocity of Snowflakes, ESA SP444 Proceedings, Millennium Conference on Antennas & Propagation, Davos, Switzerland, 9–14 April 2000.
Schönhuber, M., Lammer, G., and Randeu, W. L.: The 2Dvideodistrometer, Chapter 1, in: Precipitation: Advances in Measurement, Estimation and Prediction, edited by: Michaelides, S. C., SpringerVerlag, Berlin Heidelberg, ISBN 9783540776543, 2008.
SkofronickJackson, G., Hudak, D., Petersen, W., Nesbitt, S. W., Chandrasekar, V., Durden, S., Gleicher K. J., Huang, G., Joe, P., Kollias, P., Reed, K. A., Schwaller, M. R., Stewart, R., Tanelli, S., Tokay, A., Wang, J. R., and Wolde, M.: Global Precipitation Measurement Cold Season Precipitation Experiment (GCPEX): For Measurement's Sake, Let It Snow, B. Am. Meteorol. Soc., 96, 1719–1741, https://doi.org/10.1175/BAMSD1300262.1, 2015.
Szyrmer, W. and Zawadzki, I.: Snow Studies. Part II: Average Relationship between Mass of Snowflakes and Their Terminal Fall Velocity, J. Atmos. Sci., 67, 3319–3335, https://doi.org/10.1175/2010JAS3390.1, 2010.
Szyrmer, W. and Zawadzki, I.: Snow Studies. Part IV: Ensemble Retrieval of Snow Microphysics from DualWavelength Vertically Pointing Radars, J. Atmos. Sci., 71, 1171–1186, https://doi.org/10.1175/JASD120286.1, 2014.
Tiira, J., Moisseev, D. N., von Lerber, A., Ori, D., Tokay, A., Bliven, L. F., and Petersen, W.: Ensemble mean density and its connection to other microphysical properties of falling snow as observed in Southern Finland, Atmos. Meas. Tech., 9, 4825–4841, https://doi.org/10.5194/amt948252016, 2016.
Vega, M. A., Chandrasekar, V., Carswell, J., Beauchamp, R. M., Schwaller, M. R., and Nguyen, C.: Salient features of the dualfrequency, dualpolarized, Doppler radar for remote sensing of precipitation, Radio Sci., 49, 1087–1105, https://doi.org/10.1002/2014RS005529, 2014.
von Lerber, A., Moisseev, D., Bliven, L. F., Petersen, W., Harri, A., and Chandrasekar, V.: Microphysical Properties of Snow and Their Link to Ze–S Relations during BAECC 2014, J. Appl. Meteor. Climatol., 56, 1561–1582, https://doi.org/10.1175/JAMCD160379.1, 2017.
Westbrook, C. D. and Sephton, E. K.: Using 3Dprinted analogues to investigate the fall speeds and orientations of complex ice particles, Geophys, Res. Lett., 44, 7994–8001, https://doi.org/10.1002/2017GL074130, 2017.
Westbrook, C. D., Ball, R. C., and Field, P. R.: Radar scattering by aggregate snowflakes, Q. J. Roy. Meteor. Soc., 132, 897–914, https://doi.org/10.1256/qj.05.82, 2006.
Westbrook, C. D., Ball, R. C., and Field, P. R.: Corrigendum: Radar scattering by aggregate snowflakes, Q. J. Roy. Meteor. Soc., 134, 547–548, https://doi.org/10.1002/qj.233, 2008.
Wolfe, J. P. and Snider, J. R.: A Relationship between Reflectivity and Snow Rate for a HighAltitude SBand Radar, J. Appl. Meteor. Climatol., 51, 1111–1128, https://doi.org/10.1175/JAMCD110112.1, 2012.