There is growing interest in combining microphysical models and polarimetric radar observations to improve our understanding of storms and precipitation. Mapping model-predicted variables into the radar observational space necessitates a forward operator, which requires assumptions that introduce uncertainties into model–observation comparisons. These include uncertainties arising from the microphysics scheme a priori assumptions of a fixed drop size distribution (DSD) functional form, whereas natural DSDs display far greater variability. To address this concern, this study presents a moment-based polarimetric radar forward operator with no fundamental restrictions on the DSD form by linking radar observables to integrated DSD moments. The forward operator is built upon a dataset of >200 million realistic DSDs from one-dimensional bin microphysical rain-shaft simulations, and surface disdrometer measurements from around the world. This allows for a robust statistical assessment of forward operator uncertainty and quantification of the relationship between polarimetric radar observables and DSD moments. Comparison of “truth” and forward-simulated vertical profiles of the polarimetric radar variables are shown for bin simulations using a variety of moment combinations. Higher-order moments (especially those optimized for use with the polarimetric radar variables: the sixth and ninth) perform better than the lower-order moments (zeroth and third) typically predicted by many bulk microphysics schemes.
There is growing interest in combining numerical models and observations to further our understanding of weather and climate. For microphysical comparisons, polarimetric Doppler radar is a particularly attractive choice of observations, owing to the fact that these data can provide several key pieces of information useful for characterizing bulk properties of precipitation, such as hydrometeor sizes, shapes, concentrations, and motion (Kumjian 2013a). Radar data resolution is also ideal because it can match or exceed the higher resolution of many mesoscale and storm-scale model outputs (Keil et al. 2003). Radar data have been useful for a variety of purposes, including for model evaluation (e.g., Hagos et al. 2014; Sinclair et al. 2016; Barnes and Houze 2016; Johnson et al. 2016), data assimilation (e.g., Tong and Xue 2005; Jung et al. 2010b; Schenkman et al. 2011; Putnam et al. 2014), and gaining insights about precipitation microphysical processes (e.g., Kumjian and Ryzhkov 2010, 2012; Dawson et al. 2014; Kumjian et al. 2014; Sulia and Kumjian 2017a,b). To make such comparisons within the variable space of observations, forward operators must be used to convert model-predicted variables into quantities observed by some radar platform.
Several polarimetric radar forward operators have been developed for use with bulk and bin microphysics models (e.g., Pfeifer et al. 2008; Jung et al. 2008; Ryzhkov et al. 2011; Andrić et al. 2013). For bulk models, without exception the approach is to use model-predicted microphysical quantities to construct a particle size distribution (PSD) for each hydrometeor type at each grid box that matches the microphysics scheme’s underlying assumptions about the PSD functional form. Most often, the PSDs [including raindrop size distributions (DSDs)] are given by gamma or normalized gamma functions (e.g., Willis 1984; Testud et al. 2001). The PSD is then discretized into a series of particle size bins, whereupon electromagnetic scattering calculations are performed (e.g., Ryzhkov et al. 2011). The hydrometeor scattering properties and PSD are then integrated to obtain the radar variables of interest at each grid box.
This approach has been successful in producing simulated fields of polarimetric radar variables that reproduce basic observed signatures, particularly in convective storms (e.g., Jung et al. 2010a; Ryzhkov et al. 2011, 2013a,b; Kumjian et al. 2012, 2014, 2015; Putnam et al. 2014, 2017; Dawson et al. 2014; Johnson et al. 2016) and winter storms (e.g., Andrić et al. 2013; Sulia and Kumjian 2017a,b). However, these model–observation comparisons using forward operators often have substantial uncertainties, the details of which usually are not assessed. One class of uncertainties includes those associated with tunable parameters used or imposed by the forward operator to characterize particle properties not explicitly predicted or diagnosed by the microphysics model; for example, ice crystal shapes are treated using fixed maximum dimension–thickness relationships in the forward operators of Ryzhkov et al. (2011) and Andrić et al. (2013). Other examples include particle fall behaviors (i.e., particle orientation can strongly affect the simulated radar variables but typically is not predicted or provided by microphysics schemes), and even the choice of electromagnetic scattering calculations employed. An example of the latter is discussed by Schrom and Kumjian (2018), who quantified large errors in simulated polarimetric radar variables when homogeneous spheroids are used to approximate branched planar crystals like dendrites and stellars in scattering calculations.
Another class of uncertainties in model–observation comparison arises only when comparing to the observations themselves, originating from structural errors associated with the accuracy of approximations explicitly made in the formulation of the microphysics schemes. For example, most bulk microphysics parameterization schemes assume a functional form for the PSD, typically one that facilitates analytic integration (like the gamma PSD mentioned above). This leads to a unique mapping between model-predicted variables (e.g., total number concentration, total mass content) and radar observational quantities. Real PSDs, however, have much greater variability (i.e., not all PSDs have shapes well defined by the simple analytic PSDs assumed in most bulk schemes), leading to a greater number of degrees of freedom than bulk schemes are able to represent. In other words, although there is a unique mapping between PSD moments and radar variables for most bulk microphysics schemes, no such relationship exists in nature1 (with the exception of a nearly unique mapping between the sixth moment of the raindrop size distribution and radar reflectivity). This necessitates a treatment that accounts for this model deficiency in order to make valid comparisons between radar variables resulting from real and simulated PSDs.
This paper circumvents the problem of imposed PSD shape by creating a moment-based forward operator: one that does not assume any PSD functional form. The moment-based forward operator developed herein is flexible and can be used with a variety of bulk microphysics schemes: it directly connects the polarimetric radar variables to integrated PSD moments, regardless of the underlying PSD functional form assumed in such schemes. For example, traditional two-moment bulk microphysics schemes predict mixing ratios for mass (proportional to the third moment) and total number (the zeroth moment). Inputs from such a scheme for the forward operator are values of the zeroth and third moments at each model grid point, with no assumptions about the underlying PSD shape. Note that the moment-based approach is necessary in order to use instrument forward operators with bulk microphysics schemes that do not assume an underlying functional PSD form (e.g., Chen and Liu 2004; Szyrmer et al. 2005; Laroche et al. 2005; Kogan and Belochitski 2012). For most current bulk microphysics schemes that do use a fixed PSD functional form, our moment-based forward operator also provides an estimate of uncertainty owing to natural PSD variability not accounted for in these schemes. Such a framework is also used to find the optimal combinations of moments that minimize uncertainty in mapping to the radar variables. This can help guide the choice of prognostic variables in bulk schemes such that they are optimized for use in conjunction with radar observations.
Our study is the first step toward this moment-based approach, using the simplest framework: rain (liquid only) microphysics. Unlike the uncertainties and complexities associated with snow crystals described above, raindrop shapes are relatively well understood (e.g., Pruppacher and Pitter 1971; Beard 1976; Beard and Chuang 1987; Brandes et al. 2005; Thurai et al. 2009), as are their electromagnetic scattering properties at weather radar wavelengths (e.g., Bringi and Chandrasekar 2001; Ryzhkov et al. 2011). Additionally, dual-polarization radar variables are known to provide information—at least qualitatively—on rain microphysical processes such as evaporation (Li and Srivastava 2001; Kumjian and Ryzhkov 2010; Xie et al. 2016), size sorting (Kumjian and Ryzhkov 2012), and collision–coalescence–breakup (Kumjian and Prat 2014). In developing this moment-based forward operator, we also present the relationships between integrated DSD moments and the polarimetric radar variables, and quantify the uncertainty associated with DSD shape in terms of the polarimetric radar variables for broader use.
Section 2 outlines the methods used in this study. Section 3 describes how to identify the optimal prognostic moments for use with polarimetric radar data. The forward operator is developed in section 4. Section 5 shows example tests using a simulated rain shaft. The paper closes with a discussion and summary of the main conclusions in section 6.
DSDs are the key to linking microphysical model output to radar data because they are used to compute the bulk physical quantities of interest predicted by the model, typically total number and mass mixing ratios, as well as the radar variables, such as equivalent reflectivity factor at horizontal polarization ZH, differential reflectivity ZDR, and specific differential phase KDP. For a review of these dual-polarization radar variables, see Kumjian (2013a,b) and Kumjian (2018) and references therein.
The first step toward developing the forward operator is to create a database of DSDs. A large population of DSDs is desired because the forward operator should be able to handle any realistic precipitation situation. DSDs from both state-of-the-art bin model simulations and ground-based disdrometers are used (described below). The simulations allow for DSDs from a wide portion of the parameter space, representative of a diverse set of precipitation regimes, whereas the disdrometer data include DSDs from several different geographic regions.
The dataset will be briefly described here; details are provided in Morrison et al. (2019). Disdrometer data from the U.S. Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) program Climate Research Facility are used (Ackerman and Stokes 2003; Mather and Voyles 2013). These include samples from model RD80 Joss–Waldvogel impact disdrometers (e.g., Joss and Waldvogel 1967) and two-dimensional video disdrometers (2DVD; e.g., Tokay et al. 2001; Kruger and Krajewski 2002). The data come from geographically diverse regions, including ARM permanent sites in the U.S. Southern Great Plains, tropical western Pacific, and eastern North Atlantic (Mather and Voyles 2013; Sisterson et al. 2016; Long et al. 2016, respectively), as well as field campaigns in the Indian Ocean (Yoneyama et al. 2013; Gottschalck et al. 2013) and Finland (see Miller et al. 2016; Petäjä et al. 2016). These data cover all months and seasons, including stratiform, convective, continental, and maritime regimes. Data quality control and filtering procedures are described in Morrison et al. (2019). After these procedures, 671 303 disdrometer DSD samples remain (a sample is a 30- or 60-s average).
The simulations used herein employ the one-dimensional spectral bin microphysical model of Prat and Barros (2007) and Prat et al. (2012), following the setup used by Kumjian and Prat (2014). DSD evolution is explicitly predicted for 1 h in a one-dimensional, 3-km-tall rain shaft with 10-m vertical grid spacing. The model is initialized with a prescribed DSD at the top of the domain. At the first time step, the raindrops (which are discretized into 40 size bins) begin to fall and the DSD then freely evolves under the influence of microphysical processes. The model considers interactions among raindrops including drop coalescence, collisional and aerodynamic breakup, and sedimentation. Overall, the model performs well, consistently able to reproduce realistic DSDs as compared to radar and disdrometers (Prat et al. 2008; Prat and Barros 2009; Kumjian and Prat 2014). One potential bias exists for very heavy rainfall (>100 mm h−1), in which an overly aggressive drop breakup formulation may result in an underestimate of median drop size (Kumjian and Prat 2014). A total of 10 742 simulations were performed, covering a wide range of initial conditions, including rainfall intensity, mean drop size, and DSD shape [details can be found in Morrison et al. (2019)]. To populate the DSD dataset, DSDs are taken at every model height and output time (every 1 min). Doing so allows us to obtain samples of transient and non-steady-state DSDs from processes such as size sorting that are not well captured in disdrometer data, but are readily observed in dual-polarization radar observations (e.g., Kumjian and Ryzhkov 2012). The bin simulations produced 184 180 279 DSDs. Thus, the combined dataset is strongly dominated by the bin simulations owing to their availability.
Many of the model-predicted physical quantities of interest are proportional to specific moments of the DSD:
where Mk is the kth moment of the DSD, integrated from the minimum drop size Dmin to maximum drop size Dmax, and N(D) is the DSD (number concentration of drops with diameters in the size range D to D + dD). For example, the zeroth moment M0 of the DSD is the raindrop total number concentration, whereas the third M3 is proportional to the total raindrop mass content. For each DSD in the dataset, the integer moments k = [0, 10] were computed using Eq. (1). Because the moment values may span several orders of magnitude, we convert them to decibels (dB) using
Note that the units depend on moment order k. The moment values will be expressed in decibels for the remainder of the paper. In the current paper, we relate the polarimetric radar variables computed from observed and simulated DSDs to their respective moments. The moments themselves display natural covariability that lends itself to scaling relationships and a general DSD normalization method discussed in further detail in Morrison et al. (2019). All calculations are performed at S band (~11-cm wavelength), assuming liquid drops at 20°C and are valid for low radar antenna elevation angles (<10°). The raindrop shapes are taken as a function of size following Brandes et al. (2005). The T-matrix method (Mishchenko 2000) is used to compute the forward and backward scattering amplitudes, from which the radar variables are calculated following Ryzhkov et al. (2011). This is the same method employed by Kumjian and Prat (2014) and numerous other studies.
As mentioned in the introduction, we have not explored the effect of other sources of uncertainty such as choice of drop shape model, liquid water temperatures, and distribution of canting angles—we have focused solely on the uncertainty associated with the mapping between model-predicted quantities (integrated DSD moments) and polarimetric radar variables that is related to natural DSD variability. Uncertainty not estimated here can be easily added in subsequent work by summation of variances, assuming no correlation between different error terms. Thurai et al. (2007) showed that, at S band, the discrepancies in ZH and KDP arising from different choices of raindrop shape models and liquid water temperature are negligible, whereas ZDR differences could be up to 0.1–0.2 dB in magnitude for a small subset of DSDs characterized by large median drop sizes. For most of the DSD parameter space considered by our forward operator, then, we expect the added uncertainty arising from these choices to be smaller than the spread in ZDR values arising owing to natural variability.
Figures 1–3 show the joint histograms of M0–M10 versus the polarimetric radar variables using all DSDs in the dataset. These figures reveal the relationships between polarimetric radar variables and moments of different order k. As expected, some moments exhibit much clearer relationships with the radar variables than others. For example, ZH is nearly perfectly described by M6, whereas the dependence on M0 is rather weak (Fig. 1). This is expected given that M6 defines the radar reflectivity factor for spherical liquid drops with diameters small compared to the radar wavelength; at S band, most drops are safely considered electromagnetically small. That is to say, the Rayleigh approximation holds for all but the largest raindrops, where minor deviations from a linear ZH–M6 relationship arise. KDP (Fig. 3) appears closely related to M4 and M5 as suggested in previous studies (e.g., Sachidananda and Zrnić 1986; Bringi and Chandrasekar 2001; Lee et al. 2004; Maki et al. 2005). This is in sharp contrast to ZDR, which has more tenuous relationships with the moments, with higher-order moments displaying only slightly stronger relationships to ZDR (Fig. 2). In part, these weak relationships are because ZDR does not depend on raindrop concentration, whereas ZH and KDP do.
These joint histograms have implications for which prognostic moments offer the greatest utility for linking model output with the polarimetric radar variables. For example, M0 has a broad distribution for all three radar variables compared to higher-order moments. This implies a wide range of M0 values can produce the same ZH, ZDR, or KDP values. As such, M0 has limited utility in informing ZH, ZDR, or KDP compared to higher-order moments. From the perspective of model validation as well as data assimilation, the best prognostic DSD moments are likely to be those most directly informed by observed radar quantities; in that context, the standard choice of M0 and M3 in bulk microphysics schemes is unfortunate, as these moments are only weakly related to radar variables. This motivates the following question: which combination of moments offer the greatest information content (i.e., least spread in the joint histograms) for the radar variables? The next section addresses this question.
3. Identifying optimal predicted moments for polarimetric radar measurements
Though bulk microphysics schemes typically predict M0 and M3, the S-band radar variables are strongly related to higher moments in part because the back- and forward-scattering cross sections are proportional to D6 for particles with diameters small compared to the radar wavelength. This leads to challenges in comparing radar observations with microphysical model output. Here we assess which pair of moments (i.e., for a two-moment scheme) minimizes the variability in ZH, ZDR, and KDP for a collection of realistic DSDs, and thus would offer the most information content on those radar variables if prognosed.
We first discretize the pair of moments Mk and Mj into 1 dB × 1 dB bins. Within each bin, the ZH, ZDR, or KDP values from the DSD dataset are collected. For example, in Fig. 4a, an arbitrary bin is selected, within which M0 values range between 30 and 31 dB and M3 values between 32 and 33 dB. Within this bin, there are ~2.2 × 105 DSDs, with the bulk of their corresponding ZH values ranging from 35 to 46 dBZ (Fig. 4b). We can quantify the spread of ZH values within each pixel by calculating the standard deviation. However, before computing the standard deviation of ZH values, the evident linear trends (in logarithmic space) must be removed. Otherwise, variability within this bin will be a result of the linear trend, as opposed to the variability about this linear trend. Thus, for each 1 dB × 1 dB bin, linear trends in ZH with M0 and M3 across this bin are removed. The standard deviation is then computed, resulting in a 2D map of detrended standard deviation spanning the range of moment values (not shown).
To objectively quantify variability in prognostic moment pairs, we define a variable ξ:
where M and N are the number of bins for discretized moments Mk and Mj, respectively; is the standard deviation of the detrended polarimetric radar variable X for the mth bin of Mk and nth bin of Mj, and P is the joint normalized probability distribution function (PDF) of moments Mk and Mj in bins m and n, respectively. Physically, ξ represents the PDF-weighted spread in a given radar variable for a given pair of moments (Mk, Mj). Note that KDP is expressed in decibels for these calculations to facilitate comparison with ZH and ZDR.
The PDF weighting ensures that contributions from rare or outlier pixels are commensurate with their occurrence. However, the PDF generated by the bin simulations and disdrometer data is arbitrary (based on availability of disdrometer data and locations, choice of bin simulation parameter space, etc.) and thus may inadvertently introduce biases if used as is. Instead, a climatology of observed rainfall rates from 5-min ground-based rain gauges (see Morrison et al. 2019) is used to subsample the DSD dataset. This provides a dataset of 2 × 105 DSDs that has approximately equal contributions from the disdrometer and bin simulations and that reflects the climatological distribution of rainfall rates in the United States as measured from ground-based gauges.
The resulting ξ maps are shown in Fig. 5. One can see that ξ is minimized for different pairs of moments for ZH, ZDR, and KDP. This is expected, given that each variable has different dependencies on the DSD. For example, given that ZH is nearly equal to M6 at S band, most combinations of M6 and another moment provide the lowest ξ values (Fig. 5a). In contrast, (M5, M9) produces the lowest ξ for ZDR (Fig. 5b). Also note the large ξ values for moment order less than or equal to M3, which reveals large variability in the radar variables for the moments traditionally predicted by bulk microphysics schemes. To identify the moment pair that minimizes variability for all three variables, , , and were normalized by their respective mean values and summed together (Fig. 5d). The moment pair that yielded the minimum variability2 and thus was determined to be the optimum moment pair for informing models with dual-polarization radar observations was found to be (M6, M9). For the remainder of the paper, we will show traditional prognosed moments (M0, M3) and the ones indicated by this analysis (M6, M9). Additionally, given the practical consideration of predicting M3 in bulk microphysics schemes (as it is proportional to total mass), (M3, M6) will be shown as well. Unlike M0, M3, and M6, M9 has no conventional physical meaning3 other than the ninth moment of the DSD.
4. The moment-based forward operator
The moment-based forward operator is built using the full (combined) dataset rather than the subsampled one. This is because we desire the forward operator to cover the maximum possible spread of moment values, even if these values are rare in nature. We take a lookup table approach to the forward operator: linear interpolation (in logarithmic moment space) of the binned (Mj, Mk) values is used as a function of the input moment values. Then, the corresponding mean values (in each 1 dB × 1 dB pixel) of ZH, ZDR, and KDP are found. For example, the two-moment version of the forward operator takes as inputs a given moment pair (Mj, Mk) from, say, output from a two-moment bulk microphysics scheme that predicts Mj and Mk. The mean value for each polarimetric radar variable in the corresponding bin is assigned. Sensitivity tests (not shown) suggested 1 dB × 1 dB moment bins were an adequate balance between attaining sufficiently high resolution in the Mj–Mk parameter space and keeping the lookup tables manageable in size for our purposes herein. Note that the forward operator may be easily updated as more DSD data become available (e.g., from ongoing and future field campaigns), or with additional bin model simulations, etc., and can be generated at higher resolutions if needed in future work. A graphical depiction of three versions of the two-moment operator is shown in Fig. 6. The three versions use the (M0, M3), (M3, M6), and (M6, M9) moment pairs, respectively, for the polarimetric radar variables ZH, ZDR, and KDP. These versions of the forward operator would be used with schemes that predict M0 and M3 (most existing two-moment bulk microphysics schemes), M3 and M6, and M6 and M9, respectively.
The (M0, M3) operator (Figs. 6a–c) is based on the moments typically predicted in double-moment bulk microphysics schemes, where M0 is the total number concentration of drops and M3 is proportional to the total mass per unit volume of the drops. This version of the forward operator may be used with many commonly used two-moment bulk microphysics schemes, with the inputs simply being the predicted M0 and M3 at each model grid point. For the same M0, we see an increase in ZH, ZDR, and KDP as M3 increases. This makes sense physically: as M3 (mass) of the drops increases for a fixed number concentration, the drops must be increasing in size. A different pattern emerges for the (M3, M6) operator (Figs. 6d–f). Because M6 is almost identically ZH at S band, there is little change in ZH for increasing M3 when M6 is held fixed. For a given M3, increasing M6 leads to larger ZH, ZDR, and KDP. The (M6, M9) operator (Figs. 6g–i) is similar to the (M3, M6) operator, though a given value of the radar variables generally is spread over fewer of the 1 dB × 1 dB moment bins owing to less natural variability in the (M6, M9) moment pair (see Morrison et al. 2019).
Recall that each 1 dB × 1 dB pixel on these maps contains numerous DSDs and thus a distribution of polarimetric radar variable values within it. A novel feature of our moment-based forward operator is that it facilitates estimating the uncertainty associated with a moment-based approximation of natural DSDs. To compute such uncertainty, the standard deviation, skewness, and kurtosis are computed using the detrended data to characterize the distributions of intrinsic ZH, ZDR, and KDP variability within each pixel. Figure 7 shows standard deviation of each radar variable distribution within each pixel, for the three forward operators shown in Fig. 6. Comparing the (M0, M3) operator in the top row with the (M3, M6) and (M6, M9) in the second and third rows, a large reduction in the standard deviation of ZH is evident, which follows naturally from the fact that (M3, M6) and (M6, M9) both utilize M6. There is also a reduction in the standard deviation of ZDR and KDP evident when moving from the top to bottom rows. Skewness magnitude (Fig. 8) and kurtosis (not shown) are also substantially higher for (M0, M3), compared with the other moment-pair choices. High skewness magnitude and kurtosis imply that the uncertainty within these regions is non-Gaussian. Such non-Gaussianity poses problems for optimal estimation and Kalman-filter-based techniques that typically assume model linearity and Gaussian error statistics, with few exceptions (e.g., Hodyss 2011; Amezcua and Leeuwen 2014; Bishop 2016).
5. Examples using simulated rain shafts
To test the effectiveness of the forward operator, example rain shafts from the 1D bin simulations are used. Radar variables are calculated at each minute directly from the bin model DSDs, which are considered “truth” for these tests. We also compute the moments from these DSDs, which serve as the inputs to the forward operator. The radar variables produced by the forward operator are compared to the truth (bin simulation) values computed directly from the DSD itself. This simulation is independent from the ones used to construct the DSD database.
Figure 9 shows vertical profiles of ZH, ZDR, and KDP for a bin simulation initialized with a normalized gamma DSD aloft with rainfall rate R = 36.7 mm h−1. The “truth” profiles are shown in blue lines, whereas the (M6, M9) forward operator profiles are in gray, with plus and minus one standard deviation shown as horizontal bars on the forward-simulated profiles every 10 grid points. Each row represents a different output time in the simulation. Each profile shows the evolution of the rain shaft as raindrops fall toward the surface. At early times, size sorting of drops (e.g., Kumjian and Ryzhkov 2012; Kumjian and Prat 2014) is evident by the rapidly increasing ZDR and decreasing ZH and KDP values at the bottom edge of the rain shaft. This provides a good test for the forward operator given the somewhat exotic DSDs compared to later times when the profiles change little in height.
The forward-operator-retrieved ZH profile nearly perfectly matches the “truth” profiles at each time, which is unsurprising given that M6 is one of the moments used to inform the forward simulator. Additionally, the standard deviation is very small at all heights, as indicated by negligibly small error bars. Thus, not only does the forward operator correctly diagnose ZH, but it also correctly suggests high confidence in the diagnosis. In contrast to ZH, ZDR, and KDP provide a more difficult challenge for the forward operator given their weaker relationships to DSD moments (cf. Figs. 2 and 3). Nonetheless, the forward operator does a satisfactory job at accurately diagnosing the evolving ZDR and KDP profiles: relative error magnitudes (defined as the difference between the “truth” and forward operator curves) generally are less than 0.5%, 5%, and 10% for ZH, ZDR, and KDP, respectively. The relatively larger error magnitudes for KDP are a result of using higher-order moments (recall that M4 and M5 are the most closely related to KDP). Additionally, the diagnosed profiles are almost always within the plus and minus one standard deviation bars. For times when the “truth” lies outside the plus and minus one standard deviation bars, the diagnosed ZDR and KDP values are still well within typical theoretical radar measurement errors of ~0.1–0.2 dB for ZDR and ~0.1–0.2° km−1 for KDP (Melnikov 2004). Thus, the (M6, M9) forward operator performs well for the evolving rain shaft.
Figure 10 compares the performances of three different versions of the forward operator: (M0, M3) (Figs. 10a–c), a commonly used pair of prognostic moments for two-moment bulk microphysics parameterization schemes, (M3, M6) (Figs. 10d–f), and (M6, M9) (Figs. 10g–i). This example uses the same initial DSD aloft as in Fig. 9, with the output time t = 10 min shown. It is clear that the (M3, M6) and (M6, M9) versions are both more accurate and have less predicted spread (plus and minus one standard deviation) than the (M0, M3) version, an expected result given the higher moment orders used. In contrast, the less accurate (M0, M3) version of the forward operator has relatively larger error bars, indicating that the forward operator correctly assesses lower confidence when it is less accurate. This is a novel feature of the moment-based forward operator presented here. That the “truth” profiles fall outside the forward operator plus and minus one standard deviation bars illustrates the low-information content of M0 and M3 for the polarimetric radar variables, and is not unexpected, given that approximately 32% of all forward-simulated values will fall outside these bounds, assuming Gaussian error statistics. Furthermore, in the case of (M0, M3), Fig. 8 suggests that the standard deviation may not well characterize errors given strong deviations from Gaussianity.
Figure 11 shows another example; this time, the simulation is initialized with a normalized gamma DSD aloft with much lower rainfall rate (~0.3 mm h−1), again one that was not included in the initial dataset. As with the previous example, we see a marked improvement of the forward operator performance going from (M0, M3) (Figs. 11a–c) to (M3, M6) (Figs. 11d–f) and again to (M6, M9) (Figs. 11g–i). Once again, the plus and minus one standard deviation bars reflect the increasing forward simulator uncertainty with decreasing accuracy, particularly evident in the (M0, M3) version of the operator. The (M3, M6) operator works well for ZH and KDP, but has slight positive bias for ZDR. However, the discrepancy is ≤0.1–0.2 dB, which is well within observation error. The superior performance of the (M6, M9) operators in both examples suggests it is robust for use in both light and heavier rainfall rates. These results show that using a bulk microphysics scheme that predicts M3 and M6 (and/or M9) instead of M0 and M3 is better for use of dual-polarization radar data as a constraint or for data assimilation. Note that including M3 as a prognostic variable is important for conserving mass in models, so not predicting it may be problematic in practice. Thus, we advocate for models to use M3 and M6 as prognostic variables for two-moment schemes, or M3, M6, and M9 as prognostic variables for three-moment schemes. For two-moment bulk microphysics schemes that predict M0 and M3, dual-polarization radar data may still be used, just with considerably larger errors and greater uncertainty in the mapping between model-predicted quantities and the observed radar quantities. Whereas our forward operator attempts to quantify this uncertainty, existing forward operators use the model-assumed DSD shape (which forces a unique mapping between the predicted variables and radar variables that does not exist in nature) and do not quantify uncertainties associated with this assumption.
In principle, the approach outlined above can be extended to any number of moments and any radar variable with an accurate DSD-based forward operator. We have tested a three-moment version of the forward operator using M0, M3, and M6, the most common prognosed moments for existing three-moment schemes (e.g., Milbrandt and Yau 2005). The results showed only minimal improvement over (M3, M6) and (M6, M9) owing to the low-information content of M0 for radar variables (cf. Figs. 1–3). Higher-order moments (e.g., M3, M6, and M9) may be more useful with polarimetric radar variables and are attractive from a microphysical modeling perspective because M3 (mass) is a prognostic variable, as discussed above. Although some microphysical process rates are strongly dependent on lower-order moments, using the three-moment combination (M3, M6, and M9) allows for diagnosing lower-order moments quite well (Morrison et al. 2019) and thus is not a significant concern.
6. Discussion and summary
A large dataset of disdrometer-estimated and bin-model-simulated DSDs was constructed to quantify the relationships between different integrated moments and the S-band polarimetric radar variables, determine the uncertainty of those radar variables for a given pair of DSD moment values, and develop a moment-based polarimetric radar forward operator. This dataset comprises 671 303 DSDs estimated from Joss–Waldvogel and 2D video disdrometers at U.S. Department of Energy sites around the world, as well as 184 180 279 DSDs simulated using a one-dimensional bin microphysical model that explicitly treats raindrop collisional processes.
The data reveal a strong relationship between the sixth moment of the DSD M6 and radar reflectivity factor at horizontal polarization ZH, as expected: for spherical liquid droplets with diameters small compared to the radar wavelength, the reflectivity factor is exactly equal to M6. The specific differential phase KDP was most closely related to M4 and M5, as reported in some previous studies (e.g., Sachidananda and Zrnić 1986; Bringi and Chandrasekar 2001; Lee et al. 2004; Maki et al. 2005). In contrast, differential reflectivity ZDR showed no strong relationship with any of the DSD moments, but tended to have slightly reduced spread for higher-order moments. Future work will explore additional observations and their relationships to DSD moments, such as mean Doppler velocity from vertically pointing radar, and lidar backscatter.
The dataset was subsampled to 2 × 105 DSDs based on a climatology of observed rainfall in the United States (Morrison et al. 2019) to determine the expected natural variability of the radar variables for a given pair of moment values. The pair of moments minimizing this variability is M6 and M9. Choosing these optimal moments is a way of recasting DSD variability such that natural variability is minimized in each (Mj, Mk) pixel. In contrast, moments predicted by most bulk microphysical parameterization schemes (M0 and M3) revealed much greater variability for the polarimetric radar variables. This implies that, when comparing rain microphysical models and polarimetric radar observations, predicting higher-order moments (as opposed to or in addition to M0 and M3) could significantly improve the information content obtained from the radar variables.
A forward operator was developed to relate integrated DSD moments to polarimetric radar variables. The operator provides the mean value of ZH, ZDR, and KDP for a given pair of moment values as inputs, as well as the uncertainty in the radar variables caused by natural DSD variability (i.e., the detrended standard deviation in ZH, ZDR, and KDP within a Mj–Mk bin) and information about the distribution of radar variable values (i.e., the skewness and kurtosis). Using one-dimensional rain shafts as a benchmark, several different versions of two-moment forward operators were tested: (M0, M3), (M3, M6), and (M6, M9). The (M6, M9) version performed well for different rain shafts of varying rainfall rate, including more exotic DSDs arising from size sorting early in the rain-shaft evolution. In contrast, the forward operators with lower moment orders performed worse. The forward operator also correctly predicted its uncertainty, with greater variability indicated for the less accurate versions. This is a novel aspect of the operator developed herein.
The optimal moments for informing on the dual-polarization radar variables are of higher order than bulk microphysics schemes typically prognose. Though such high moments individually may not provide much of a constraint for lower-order moments needed for such schemes, they can still reduce the uncertainty considerably when used in combination (i.e., multiple prognostic moments) and/or in combination with lower-order moments (Morrison et al. 2019). In other words, if attempting to diagnose the kth moment Mk, reference moment Mk+n always provides a better estimate than reference moment Mk−n for all n. Further, Morrison et al. (2019) show that a three-moment normalization using M3, M6, and M9 will result in only ~21% of the variability in M0 compared to not using the DSD normalization. Thus, use of such higher-order moments in bulk microphysics schemes may not be detrimental, and indeed could be beneficial when combined with lower-order moments typically prognosed (like M3).
The moment-based forward operator developed herein is necessary for coupling radar observations with bulk microphysics schemes that do not assume a DSD functional form (e.g., Chen and Liu 2004; Szyrmer et al. 2005; Laroche et al. 2005; Kogan and Belochitski 2012). Other forward operators reliant on a discretized DSD would require assuming an explicit DSD functional form, imposing structural error into the mapping between model output and radar observations. The approach herein strives to minimize and quantify this type of uncertainty, such that the majority of the remaining uncertainty contained in the forward operator arises owing to DSD natural variability, and is explicitly estimated. Ultimately, this type of approach should lead to improved mapping of model output to the radar observational parameter space with a better characterization of uncertainty. For traditional bulk microphysics schemes, use of the moment-based operator described here may prevent errors associated with overconfidently comparing approximate bulk schemes to observations associated with more complex, realistic DSDs.
This work was funded by U.S. DOE Atmospheric System Research Grant DE-SC0016579. The National Center for Atmospheric Research is sponsored by the National Science Foundation. Disdrometer data were obtained from the Atmospheric Radiation Measurement (ARM) Climate Research Facility Data Archive in Oak Ridge, Tennessee, United States, compiled and maintained by M. Bartholomew. Argonne National Laboratory’s work was supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, under Contract DE-AC02-06CH11357. We thank the anonymous reviewers for their helpful suggestions and comments that improved the clarity of the manuscript.
Note that these uncertainties are also relevant to the calculation of microphysical process rates; however, this is beyond the scope of the current study and will be addressed in future work.
Note that this is somewhat sensitive to how the ξ for each variable are summed. Different weightings may be applied as needed. For example, if less confidence is placed on comparing ZDR to observations owing to calibration issues, or on KDP owing to difficulties in its estimation because of noisy total differential phase fields, one could weight the summation away from one of the variables in favor of the other two.
If normalized by M6, then M9 could be considered the “reflectivity-weighted mass” of the distribution.