Abstract

In radar polarimetry, the differential phase ΨDP consists of the propagation differential phase ΦDP and the backscatter differential phase δhv. While ΦDP is commonly used for attenuation correction (i.e., estimation of the specific attenuation A and specific differential phase KDP), recent studies have demonstrated that δhv can provide information concerning the dominant size of raindrops. However, the estimation of ΦDP and δhv is not straightforward given their coupled nature and the noisy behavior of ΨDP, especially over short paths. In this work, the impacts of estimating ΦDP on the estimation of A over short paths, using the extended version of the ZPHI method, are examined. Special attention is given to the optimization of the parameter α that connects KDP and A. In addition, an improved technique is proposed to compute δhv from ΨDP and ΦDP in rain. For these purposes, diverse storm events observed by a polarimetric X-band radar in the Netherlands are used. Statistical analysis based on the minimum errors associated with the optimization of α and the consistency between KDP and A showed that more accurate and stable α and A are obtained if ΦDP is estimated at range resolution, which is not possible by conventional range filtering techniques. Accurate δhv estimates were able to depict the spatial variability of dominant raindrop size in the observed storms. By following the presented study, the ZPHI method and its variations can be employed without the need for considering long paths, leading to localized and accurate estimation of A and δhv.

1. Introduction

Conventional S- and C-band weather radars have been used for several decades to monitor the evolution of precipitation. In recent years the technology of those conventional radars has been upgraded to polarimetric technology in order to further improve weather radar measurements (Doviak et al. 2000). Severe weather can produce rapid and localized surface damage associated with, for example, heavy rain and tornadoes. In this context, a network of small polarimetric X-band weather radars may be suitable to obtain observations of fast-developing storms at close range and at resolutions higher than those from conventional radars (McLaughlin et al. 2009; Chandrasekar et al. 2018).

One of the advantages of polarimetric radars is given by the measurements of differential phase between the horizontally and vertically polarized signals caused by the delay of one with respect to the other as both signals propagate through hydrometeors. In this way, the differential phase ΨDP (°) is independent of attenuation, miscalibration, and partial beam blockage (PBB) effects (Doviak and Zrnić 1993). However, ΨDP measurements can include phase shifts in the backward direction as a result of Mie scattering, the so-called backscatter differential phase δhv (°), and random fluctuations ε (°) on the order of few degrees. In general, a ΨDP range profile is modeled as

 
ΨDP(r)=ΦDP(r)+δhv(r)+ε,
formula

where ΦDP(r) (°) represents the differential phase in the forward direction and r (km) indicates the distance from the radar. Two useful variables that can be estimated from ΦDP are the specific differential phase KDP (° km−1) and the specific attenuation A (dB km−1), which are commonly used for the estimation of rainfall rate and attenuation correction (Bringi and Chandrasekar 2001).

The traditional method to estimate KDP (or ΦDP) from ΨDP when δhv is significant is given by Hubbert and Bringi (1995), and several attempts have been proposed to improve KDP estimates at X-band frequencies (Wang and Chandrasekar 2009; Giangrande et al. 2013; Schneebeli et al. 2014; Huang et al. 2017). The specific differential phase KDP has been used to correct measurements of reflectivity Z (dBZ) affected by radar calibration and PBB (Giangrande and Ryzhkov 2005). In addition, KDP has led to improved estimation of rainfall rate, mostly in heavy rain or mix rain, because of its quasi-linear relation to liquid water content (Lim et al. 2013). Although radar measurements seem to benefit from using KDP, comprehensive research on KDP is still needed because it is a challenge to provide accurate KDP from noisy measurements of ΨDP.

Existing methods to estimate A in rain assume that A = αKDP, where α is a constant for a given frequency (Bringi et al. 1990). Testud et al. (2000) also used the relation between A and KDP in their rain profiling ZPHI technique, to express A in terms of the difference of ΦDP and measurements of Z, avoiding KDP calculation. However, it is known that α is sensitive to temperature, drop size distribution (DSD), and drop size variabilities; therefore, Bringi et al. (2001) extended the ZPHI technique to avoid a priori value for α. These methods have been adapted to address attenuation problems at X-band frequencies (Matrosov et al. 2002; Park et al. 2005a; Gorgucci et al. 2006; Lim and Chandrasekar 2016). Moreover, Ryzhkov et al. (2014), Wang et al. (2014), and Diederich et al. (2015) modified the extended ZPHI method to improve rainfall-rate estimation and to demonstrate that A can be used to reduce issues related to radar calibration and PBB. Despite these promising benefits, the potential of using A might be limited depending on the approach to obtain ΦDP and α (Bringi et al. 2001; Ryzhkov and Zrnić 2005).

In contrast to KDP and A, limited research has been conducted on the applications of δhv. For example, δhv can be a suitable candidate to mitigate uncertainties related to the differential reflectivity ZDR (dB) because δhv and ZDR offer a correlated behavior (Scarchilli et al. 1993; Testud et al. 2000) and because δhv is independent of attenuation and radar calibration; see Eq. (1). These aspects of δhv could be useful to establish relations between δhv and the median drop diameter D0 (mm) (Trömel et al. 2013) because D0 is often expressed in terms of ZDR (Matrosov et al. 2005; Kim et al. 2010). Moreover, Otto and Russchenberg (2010) included δhv estimates to retrieve DSD parameters. Hubbert and Bringi (1995), Otto and Russchenberg (2011), and Trömel et al. (2013) estimated δhv by subtracting ΦDP from ΨDP, while Schneebeli and Berne (2012) included a Kalman filter approach. The effectiveness of estimating δhv at high resolution is rather complicated because of the cumulative and noisy nature of ΨDP and possible remaining fluctuations on ΦDP.

The purpose of this work is to 1) explore the role and impact of estimated ΦDP profiles on the performance of the extended ZPHI method at X-band frequencies to improve estimates of α and A over short paths and 2) develop a technique to compute δhv in rain while keeping the spatial variability of drop sizes. For such purpose, two KDP (or ΦDP) methods, by Hubbert and Bringi (1995) and Reinoso-Rondinel et al. (2018), are reviewed in section 2 as well as three attenuation correction approaches, by Bringi et al. (1990), Testud et al. (2000), and Bringi et al. (2001). In addition, the δhv algorithm is introduced, which integrates estimates of KDP and A. In section 3, the performances of the attenuation correction methods that assume a constant α are compared using four storm events. This comparison is extended in section 4 to examine the selection of α profile by profile and its impact on A and Z. In section 5, the δhv technique is evaluated. Section 6 focuses on the statistics of α, A, Z, and δhv to conduct further assessments of the presented methods. Finally, section 7 draws conclusions of this article.

2. Estimation techniques for ΨDP-based variables

a. Estimation of KDP

In the conventional technique given by Hubbert and Bringi (1995), a low-pass filter is designed such that gate-to-gate fluctuations at scales of the range resolution Δr (km) are filtered from a ΨDP(r) profile. Fluctuations at range scales larger than Δr (i.e., δhv “bumps”) are removed by applying the same filter multiple times to new generated ΨDPg profiles by combining a previous filtered and original ΨDP profile. In this manner the corresponding ΦDP profile is obtained and KDP is given by taking a range derivative of ΦDP. For the generation of ΨDPg, a predetermined threshold τ (°) is required, which is on the order of 1–2 times the standard deviation of ΨDP, hereafter σP (°). One of the limitations of this technique is that accurate estimates of ΦDP and KDP at Δr scales are hardly achieved (Grazioli et al. 2014).

An adaptive approach that estimates KDP at high spatial resolution while controlling its standard deviation σK (° km−1) is given by Reinoso-Rondinel et al. (2018). For notation purposes, the difference of a radar variable V over a given pathlength is expressed as ΔV. Besides ΨDP, attenuation-corrected Z and ZDR profiles are also required, as well as a predefined pathlength interval [Lmin;Lmax] (km). For gate i, located at range ri, a set of σK samples are obtained from [Lmin;Lmax] using a theoretical expression of σK. The pathlength that minimizes the σK set is selected and denoted as L(i). Assuming the correlated behavior between ZDR and δhv, ΔΨDP samples in the range [riL(i); ri+L(i)] that do not satisfy the condition |ΔZDR|<σZDR are filtered to avoid contamination from Δδhv. The standard deviation of the ZDR profile is denoted as σZDR. The spatial variability of ΨDP at Δr scales is captured by downscaling each remaining ΔΨDP sample from L(i) to Δr scale. A downscaling parameter w(i) [0, 1] is derived from Z and ZDR in the same interval [riL(i); ri+L(i)], and KDP(i) is estimated as

 
KDP(i)=1Mj=1MΔΨDP(j)w(j)(i)2Δr,withj=1,2,,M,
formula

where M represents the number of ΔΨDP samples with negligible Δδhv. The actual σK(i) is calculated using the terms inside the sum operation in Eq. (2). The KDP and σK profiles are obtained by repeating the same procedure over the remaining gates, while the corresponding ΦDP profile is calculated by simply integrating KDP in range. In addition, a profile of the normalized standard error (NSE) of KDP is given by the ratio between actual σK and KDP. This approach was demonstrated for rain particles at X-band frequencies, and therefore any undetected Z and ZDR echoes from hydrometeors other than rain can lead to inaccurate KDP estimates. The two KDP methods will be referred to as the conventional (C) and the adaptive high-resolution (AHR) approaches, respectively. A diagram is presented in Fig. 1 to briefly indicate the inputs and outputs of each method.

Fig. 1.

Methods associated with the estimation of (a) KDP and (b) A. The outputs related to the conventional KDP technique are indicated with red, while the outputs related to the adaptive high-resolution approach are indicated with green.

Fig. 1.

Methods associated with the estimation of (a) KDP and (b) A. The outputs related to the conventional KDP technique are indicated with red, while the outputs related to the adaptive high-resolution approach are indicated with green.

b. Estimation of A

For attenuation correction purposes, Z and ZDR profiles are represented as Z(r)=z(r)+PIA(r) and ZDR(r)=zdr(r)+PIADP(r), respectively, where z (dBZ) and zdr (dB) represent the attenuated reflectivity and the attenuated differential reflectivity, respectively; and PIA(r) (dB) indicates the two-way path-integrated attenuation in reflectivity and PIA(r)DP (dB) in differential reflectivity.

Bringi et al. (1990) introduced the differential phase (DP) approach such that A(r) = αKDP(r) and PIA(r)=αΦDP(r), where α [dB (°)−1] is assumed to be a constant coefficient. Gorgucci and Chandrasekar (2005) studied the accuracy of this method using simulated radar variables at X-band frequencies and showed that estimates of A are very sensitive to inaccurate estimates of KDP, while estimates of PIA lead to Z values associated with only a slight degradation of the average error for attenuation correction, ±1.5 dB.

To improve the DP method, Testud et al. (2000) introduced the ZPHI method that estimates A(r) in a path interval [rp;rq], where rq>rp. First, A(r) is expressed as a function of two known variables, z(r) and z(rq), and one unknown, A(rq). Then, A(rq) is obtained using z(rq) and the empirical relation ΔPIA = αΔΦDP, where ΔPIA = PIA(rq) − PIA(rp) and ΔΦDP = ΦDP(rq)ΦDP(rp). In this way, A(r) is estimated at Δr scales, reducing errors related to KDP(r). Although [rp;rq] can be freely selected; ΔΦDP could be inaccurate at short path intervals and/or be contaminated by δhv(rp) and δhv(rq). In addition, if z(r) includes localized observations of hail or mixtures of rain and hail in [rp;rq], then A(r) might be biased over the entire path interval.

Using a constant α may lead to limited approximations of A(r) and PIA(r) because α is sensitive to DSD, drop shape, and temperature variabilities (Jameson 1992). To take into account the sensitivity of α, Bringi et al. (2001) extended the ZPHI method to search for optimal α values at C-band frequencies, called the CZPHI method. An initial value for α is selected from a predefined interval [αmin;αmax], and A(r) is estimated according to the ZPHI method. The estimated A(r) is integrated over [rp;rq] to build a differential phase profile denoted as ΦDP(r,α). Repeating this procedure for the remaining values of α, the optimal α is the one that minimizes the error E (°) given by

 
E=i=pq|ΦDP(ri,α)ΦDP(ri)|,withi=p,,q.
formula

Note that the optimization process requires the estimation of ΦDP, which implies the need for a proper way to filter noise and δhv components from ΨDP while maintaining its spatial variability. However, meeting such requirements is not straightforward; therefore, the reliability of an “optimal” α to estimate A and PIA depends on the performance of the chosen approach to estimate ΦDP. The inputs and outputs associated with the three presented attenuation correction methods are summarized in Fig. 1.

To determine PIA(r)DP, integrate the specific differential attenuation ADP(r) (dB km−1) that is given by ADP=γA. The DP and ZPHI methods assume γ to be constant, whereas the CZPHI technique searches for an optimal γ, addressing its sensitivity to DSD variability (i.e., rain type). However, such sensitivity of γ is less at X-band frequencies than at C- and S-band frequencies (Ryzhkov et al. 2014). In this work, ADP will be given by ADP = γA(CZPHI), where A(CZPHI) represents the specific attenuation determined by the CZPHI approach and γ is assumed a constant.

Representative values for α and γ at X-band frequencies can be given by the mean fit of simulated polarimetric relations using a large set of DSDs and different drop shapes and temperatures. For example, Kim et al. (2010) and Ryzhkov et al. (2014) demonstrated that α values vary in the interval [0.1; 0.6] dB (°)−1, and Otto and Russchenberg (2011) obtained an average value of 0.34 dB (°)−1 for α and for γ a value of 0.1618. Similar results were suggested by Testud et al. (2000), α=0.315 dB (°)−1; Kim et al. (2010), α=0.35 dB (°)−1; and Snyder et al. (2010) ,α=0.313 dB (°)−1; while Ryzhkov et al. (2014) estimated γ equal to 0.14 for tropical rain (i.e., low ZDR and high KDP) and 0.19 for continental rain (i.e., high ZDR and low KDP). It is important to note that other authors have suggested smaller average values for α. For example, Bringi and Chandrasekar (2001) simulated polarimetric variables in rain and indicated that α=0.23 dB (°)−1. Matrosov et al. (2014) avoided simulations by using observations resulting from collocated X- and S-band radars and found α in the range of 0.20–0.31 dB (°)−1. Thus, a representative value for α can vary depending on models and assumptions used to simulate polarimetric variables, on the type of observed storms and their geographical locations, and on the accuracy of measurements.

c. Estimation technique for δhv

A δhv approach is presented to identify and separate Mie scattering signatures from noise and random fluctuations embedded in ΨDP. A flowchart of the δhv algorithm is illustrated in Fig. 2. Three inputs are required: a 2D ΨDP field measured in rain, the corresponding KDP field obtained from the AHR approach, and the A field estimated by the CZPHI method. Given these inputs, the resulting δhv field is based on the following five steps:

  1. Design and apply a filter to smooth strong outliers from a ΨDP profile, taking Δr into account. Correct each smoothed ΨDP profile for system phase offset by subtracting the mean of ΨDP over the first 5% of measured gates.

  2. Obtain ΦDP by integrating profiles of A, if they are associated with a minimum error E, otherwise by integrating KDP profiles. Next, subtract ΦDP from ΨDP, profile by profile, as a first attempt to estimate the corresponding δhv field. The next steps are related to 2D processing.

  3. Remove unusual δhv values larger than 12° from the δhv field. According to Testud et al. (2000), Trömel et al. (2013), and Schneebeli et al. (2014), the simulated δhv values at X-band frequencies rarely reach 12°. The remaining noise in δhv is reduced by assuming that similar values of δhv are collocated with similar values of KDP as follows. Set Kmin as the minimum of KDP and Kmax as Kmin+ΔK, where ΔK (° km−1) is given by Eq. (4). Define S as a set of δhv samples, whose gates are collocated with KDP values in the interval [Kmin;Kmax]. Reject δhv samples from S that are outside the interval [δ¯hvυσδhv; δ¯hv+υσδhv], where δ¯hv and σδhv indicate the arithmetic mean and the standard deviation of the samples in S, respectively; υ is a predefined threshold in the interval [1;2] and a value of 1 is chosen. This process is iterated by shifting [Kmin;Kmax] toward high values in small steps such that Kmin = Kmax and Kmax = Kmin+ΔK until Kmax is equal to the maximum of KDP. To obtain sufficient samples in S, ΔK is given as 
    ΔK={0.2Kmin2.5°km1,0.52.5°<Kmin<8°km1,1.0Kmin8°km1,
    formula
    because high KDP values are less frequent than small KDP values (e.g., see the KDP fields in Figs. 3, 8, and 11).
  4. Apply a 2D interpolation method to fill empty gaps on δhv caused by step 3. For this task, the inpainting (or image fill-in) algorithm (Bertalmio et al. 2003; Criminisi et al. 2004; Elad et al. 2005) is selected because it is one of the image processing algorithms commonly used to smoothly interpolate 2D images. The essential idea is to formulate a partial differential equation (PDE) for the “hole” (interior unknowns) and to use the perimeter of the hole to obtain boundary values. The solution for the interior unknowns involves the discretization of PDEs on the unknowns’ points into a system of linear equations. D’Errico (2006) implemented an inpainting code for 2D arrays that is freely available and used for this step. The code offers multiple methods to formulate a PDE, and the method referred to as the spring method is selected because it provides a reasonable compromise between accuracy and computational time.

  5. (optional) To better distinguish storm cells from their background (i.e., for radar displaying purposes), it is recommended to replace areas of δhv that are linked to |KDP| < 0.4° km−1 (i.e., weak rain echoes) by a representative value. This value is chosen as the mean of δhv samples constrained by |KDP| < 0.4° km−1 and |δhv|<σ¯δhv, where σ¯δhv indicates the mean of σδhv samples obtained in a similar manner as in step 3 but using δhv after step 4. The value of 0.4° km−1 is found to match the 30-dBZ level used in this work for storm cell identification.

Fig. 2.

A flowchart for the estimation of δhv. It consists of five steps, where steps 1 and 2 are processed in 1D (i.e., along a PPI radial), while steps 3–5 are processed in 2D (i.e., a complete PPI).

Fig. 2.

A flowchart for the estimation of δhv. It consists of five steps, where steps 1 and 2 are processed in 1D (i.e., along a PPI radial), while steps 3–5 are processed in 2D (i.e., a complete PPI).

Fig. 3.

Observations by IDRA radar at elevation angle of 0.5° in the NL at 1216 UTC 18 Jun 2011, event E1. Fields of (a) differential phase ΨDP, (b) z, (c) KDP(C) from the conventional approach, and (d) KDP(AHR) from the AHR approach. In (b)–(d), attenuation-corrected 30-dBZ levels are indicated by black contour lines; in (c) and (d), −1° km−1 levels are indicated by magenta contour lines. The red rings are at 5-km increments.

Fig. 3.

Observations by IDRA radar at elevation angle of 0.5° in the NL at 1216 UTC 18 Jun 2011, event E1. Fields of (a) differential phase ΨDP, (b) z, (c) KDP(C) from the conventional approach, and (d) KDP(AHR) from the AHR approach. In (b)–(d), attenuation-corrected 30-dBZ levels are indicated by black contour lines; in (c) and (d), −1° km−1 levels are indicated by magenta contour lines. The red rings are at 5-km increments.

3. Evaluation of KDP processing by the ZPHI method

a. Datasettings and preprocessing

The polarimetric X-band International Research Center for Telecommunications and Radar (IRCTR) Drizzle Radar (IDRA; Figueras i Ventura 2009) is located at the Cabauw Experimental Site for Atmospheric Research (CESAR) observatory in the Netherlands (NL) at a height of 213 m from ground level (Leijnse et al. 2010). Its operational range and range resolution are equal to 15.3 and 0.03 km, respectively, while the antenna rotates over 360° in 1 min. Four storm events, E1–E4, that occurred in the Netherlands during the year 2011 will be used for demonstration and analysis purposes. A description of these events is summarized in Table 1.

Table 1.

Description of four storm events E1–E4 observed in the Netherlands.

Description of four storm events E1–E4 observed in the Netherlands.
Description of four storm events E1–E4 observed in the Netherlands.

To remove areas that include particles other than rain and/or areas with low signal-to-noise ratio (SNR), measurements of linear depolarization ratio LDR (dB) are used, such that range gates with LDR larger than −18 dB are discarded from ΨDP, z, and zdr fields. Further preprocessing includes suppressing isolated segments of a ΨDP profile smaller than 0.25 km and rejecting a ΨDP profile if the percentage of gates with measurements is less than 5%. Because a ΨDP profile could be noisy at ranges behind strong reflectivity echoes associated with low SNR and fully attenuated signals, its range extent needs to be determined. The ending range of a ΨDP profile is determined based on σ¯P, which represents the average of multiple σP samples by running a five-gate window along the ΨDP profile. If σ¯P is less than 1.5°, then the ending range is given by the last measured gate in the downrange direction. Otherwise, the ending range is set by the middle gate of the second consecutive window whose σP values are less than σ¯P, starting at the last measured gate and moving toward the radar. The ending range is used to limit the corresponding extent of z and zdr profiles. After this, σ¯P is calculated again to estimate KDP by the conventional technique.

b. Comparison between KDP and A

Next, KDP(C) and KDP(AHR) will be compared against A(ZPHI) using the empirical relation A = αKDP, where α is 0.34 dB (°)−1, as suggested by Otto and Russchenberg (2011). In this scheme, A(ZPHI) is used as a reference to evaluate both KDP techniques and their impact on Z.

To estimate KDP(C), a finite impulse response (FIR) filter is used such that the order of the filter is 36 and the cutoff range scale is 1 km, including a Hann window. The required threshold τ is set to 1.5σ¯P. Such a filter design is found suitable for Δr= 0.03 km. For the estimation of KDP(AHR), values of L on the order of 3 km are associated with theoretical values of σK< 0.5° km−1 for Δr = 0.03 km (Reinoso-Rondinel et al. 2018) and therefore [Lmin;Lmax] is predefined as [2;5] km. The z and zdr inputs are corrected for attenuation and differential attenuation, respectively, according to the DP method, in which a linear regression fit of 1 km is applied to ΨDP profiles. To estimate σZDR a five-gate window is run along a given ZDR profile. For the calculation of A(ZPHI), ΔΦDP is derived from ΦDP(C) instead of ΦDP(AHR) to evaluate KDP(AHR) in an independent manner. A path interval [rp;rq] is defined by the first and last data points, in the downrange direction, of a ΦDP(C) profile. In cases where ΔΦDP< 0° as a result of a reduced SNR profile, the estimation of A(ZPHI) is avoided.

Results from the storm event E1 at 1216 UTC are shown in Fig. 3. The ΨDP field shows a rapid increment in range on the north side of the storm, whereas ΨDP rarely increases on the south side. Note that the ΨDP field is not adjusted for phase offset. The attenuated z field represents a relatively small cell of a nonuniform structure in close proximity to the radar. The 30-dBZ contour is obtained from the attenuation-corrected Z using the ZPHI method [i.e., after calculating A(ZPHI) as explained previously]. Comparing KDP(C) and KDP(AHR), the KDP(AHR) field is able to maintain the spatial variability of the storm down to range resolution scale, eliminating areas of KDP smaller than −1° km−1, which are present in KDP(C). However, the coverage of the KDP(AHR) field is smaller than that of KDP(C). This is because in the AHR approach, it is not always possible to obtain ΔΨDP samples with negligible Δδhv; that is, M=0 in Eq. (2). Note that isolated KDP segments smaller than 2 km were removed from both KDP fields in order to avoid estimates of KDP that could be associated with noisy areas and/or low accuracy.

The scatterplots KDP(C)–A(ZPHI) and KDP(AHR)–A(ZPHI) resulting from the same event, E1, are compared in Fig. 4. In Fig. 4a, it can be seen that the KDP(AHR)–A(ZPHI) scatterplot (14 783 data points) is more consistent than that of KDP(C)–A(ZPHI) (15 490 data points) with respect to the empirical relation A = 0.34KDP. In a quantified comparison, the correlation coefficient ρKA between KDP(C) and A(ZPHI) is equal to 0.65, whereas for KDP(AHR) and A(ZPHI) it is 0.96. Their corresponding standard deviations σKA with respect to the empirical relation are 1.20 and 0.41° km−1, respectively. To compare the impact of both KDP techniques on the DP method, z values are corrected for attenuation using the DP and ZPHI correction methods, and are denoted as Z(DP, C), Z(DP, AHR), and Z(ZPHI, C); see Fig. 1. The scatterplots Z(DP, C)–Z(ZPHI, C) and Z(DP, AHR)–Z(ZPHI, C) are compared in Fig. 4b such that Z(ZPHI, C) estimates are used as reference. It is observed that for relatively high values of Z(ZPHI, C), Z(DP, C) values are slightly overcorrected, which agrees with Gorgucci and Chandrasekar (2005) and Snyder et al. (2010). In contrast, Z(DP,AHR) values are found significantly consistent with Z(ZPHI, C) estimates. The mean biases associated with Z(DP,C) and Z(DP, AHR) are equal to 0.95 and −0.21 dB, respectively, for Z(ZPHI, C) ≥35 dBZ. The errors quantified by ρKA, σKA, and bias Z are summarized in Table 2. The remaining events, E2–E4, at 1450, 1955, and 0558 UTC, respectively, were also analyzed in a similar manner and the corresponding quantified errors are indicated in Table 2.

Fig. 4.

(a) The KDP(C)–A(ZPHI) scatterplot resulting from event E1 at 1216 UTC is indicated by red dots, and the KDP(AHR)–A(ZPHI) scatterplot is indicated by green dots. In addition, the empirical relation KDP = (1/α)A is indicated by the black line, where α = 0.34 dB (°)−1. (b) As in (a), but for Z(DP, C)–Z(ZPHI) and Z(DP, AHR)–Z(ZPHI) scatterplots. Also, the relation Z(DP) = Z(ZPHI) is indicated by the black line. The biases are computed for Z(ZPHI) ≥ 35 dBZ.

Fig. 4.

(a) The KDP(C)–A(ZPHI) scatterplot resulting from event E1 at 1216 UTC is indicated by red dots, and the KDP(AHR)–A(ZPHI) scatterplot is indicated by green dots. In addition, the empirical relation KDP = (1/α)A is indicated by the black line, where α = 0.34 dB (°)−1. (b) As in (a), but for Z(DP, C)–Z(ZPHI) and Z(DP, AHR)–Z(ZPHI) scatterplots. Also, the relation Z(DP) = Z(ZPHI) is indicated by the black line. The biases are computed for Z(ZPHI) ≥ 35 dBZ.

Table 2.

Comparison results between KDP(C) estimates and KDP(AHR) using as a reference values of A(ZPHI) resulting from the ZPHI method for four storm events. Data points in each event are given: event E1 (~14 000), E2 (~13 000), E3 (~40 000), and E4 (~30 000).

Comparison results between KDP(C) estimates and KDP(AHR) using as a reference values of A(ZPHI) resulting from the ZPHI method for four storm events. Data points in each event are given: event E1 (~14 000), E2 (~13 000), E3 (~40 000), and E4 (~30 000).
Comparison results between KDP(C) estimates and KDP(AHR) using as a reference values of A(ZPHI) resulting from the ZPHI method for four storm events. Data points in each event are given: event E1 (~14 000), E2 (~13 000), E3 (~40 000), and E4 (~30 000).

From the previous analysis, the following can be highlighted. The values of KDP(AHR) and A(ZPHI), determined by two independent methods, show a strong agreement to the empirical relation A=αKDP, leading to equivalent Z(DP,AHR) and Z(ZPHI,C) results. In the contrary, the agreement between KDP(C) and A(ZPHI) is less evident, and although KDP(C) barely includes substantial errors on attenuation-corrected Z(DP,C), it can significantly impact estimates of A by the DP method. Similar findings at X-band frequencies were reported by Gorgucci and Chandrasekar (2005) but using simulated data.

4. Impact of KDP processing on the CZPHI method

In this section, the ability to estimate ΦDP by both KDP approaches is studied and their impact on the performance of finding optimal α values for the estimation of A and the correction of Z by the CZPHI method is measured. For analysis purposes, the minimum E obtained from Eq. (3) is expressed as E = ei, with i = p,,q, where ei represents the minimum error at range ri. As such, the arithmetic mean and standard deviation of ei, e¯min (°) and σemin (°), respectively, will be used as quality measures.

At X-band frequencies, [αmin;αmax] is predefined as [0.1;0.6] dB (°)−1 with steps of 0.02 dB (°)−1, as suggested by Park et al. (2005b) and Ryzhkov et al. (2014). For a correct optimization process, it is recommended that rqrp should be at least 3 km and that ΔΦDP be larger than 10°. In addition, if the ΦDP(C) profile is used in Eq. (3), then the percentage of gates with KDP> 0° km−1 should be at least 50%, whereas if the ΦDP(AHR) profile is used, the percentage of gates with KDP> 0.5° km−1 and NSE < 20% should be larger than 80%. The percentage threshold for ΦDP(C) is less than for ΦDP(AHR) because the conventional method rarely avoids negative KDP values. If these conditions are met, α is selected by minimizing E, considering only range gates that satisfy the stated conditions; otherwise α is equal to 0.34 dB (°)−1.

a. Event E1: Single cell

1) Optimization analysis

Results involved in the optimization process along azimuth 288.1° for storm event E1 at 1216 UTC are shown in Figs. 5a–c. In Fig. 5a, it is seen that the minimum E when ΦDP(C) is used is much larger than when ΦDP(AHR) is used and their corresponding optimal values for α are αΦDP(C) = 0.24 and αΦDP(AHR) = 0.34 dB (°)−1. The reason why the two α values are different can be explained by observing the measured ΨDP and the estimated ΦDP(C) and ΦDP(AHR) profiles shown in Figs. 5b and 5c, respectively. First, note that ΨDP might include (i) a δhv bump in the range [3.5;5.5] km and (ii) oscillations in the range [6.5;8.5] km. Second, the δhv bump is more noticeable in ΦDP(C) than in ΦDP(AHR). In consequence, the matching between ΦDP(C) and ΦDP(CZPHI) shown in Fig. 5b is not as good as the one observed in Fig. 5c. Note that ΦDP(CZPHI) represents ΦDP(ri,α) in Eq. (3). The extent of the ΦDP(AHR) profile is less than that of ΦDP(C) because M in Eq. (2) appears to be 0 at the beginning and ending ranges of ΨDP. However, this limited extent of ΦDP(AHR) avoids the oscillations seen at the ending ranges of ΨDP.

Fig. 5.

(a) Errors obtained from Eq. (3): EΦDP(C) (red) and EΦDP(AHR) (green) in the azimuth 288.1°, as a function of α [αmin;αmax]. (b) Profiles of ΨDP, ΦDP(C), and ΦDP(CZPHI) are shown as a function of range. In addition, upper and lower ΦDP(CZPHI) bounds (dashed lines) corresponding to αmin and αmax, respectively. (c) As in (b), but using ΦDP(AHR) rather than ΦDP(C). (d) Stemplots of selected αΦDP(C) (red) and αΦDP(AHR) (green) as a function of azimuth.

Fig. 5.

(a) Errors obtained from Eq. (3): EΦDP(C) (red) and EΦDP(AHR) (green) in the azimuth 288.1°, as a function of α [αmin;αmax]. (b) Profiles of ΨDP, ΦDP(C), and ΦDP(CZPHI) are shown as a function of range. In addition, upper and lower ΦDP(CZPHI) bounds (dashed lines) corresponding to αmin and αmax, respectively. (c) As in (b), but using ΦDP(AHR) rather than ΦDP(C). (d) Stemplots of selected αΦDP(C) (red) and αΦDP(AHR) (green) as a function of azimuth.

The selected αΦDP(C) and αΦDP(AHR) values as a function of azimuth for the same storm are depicted in Fig. 5d. Values for α that are related to a minimum E (i.e., optimal α values) are encircled by black edges, while those that are nonrelated to a minimum E are represented without edges. Note that optimal αΦDP(AHR) values are close to 0.34 dB (°)−1, whereas those related to ΦDP(C) are mostly smaller than 0.34 dB (°)−1 and sometime equal to αmin. An optimal α that equals αmin or αmax could be associated with an inadequate matching between the input ΦDP and the obtained ΦDP(CZPHI), which can lead to incorrect α. The resulting e¯min values associated with ΦDP(C) and ΦDP(AHR) are 2.16° and 0.20°, respectively, and their corresponding σemin values are 0.75° and 0.08°. These results come from the azimuthal sector [280°;310°], which covers approximately the north side of the storm shown in Fig. 3. Outside this sector, the constant α was selected, associated with either ΦDP(C) or ΦDP(AHR), because the stated conditions were not met.

2) Performance analysis

The impact of the optimal selection of αΦDP(C) and αΦDP(AHR) on the estimation of A(CZPHI) is measured using KDP(AHR) as a reference because of 1) the consistency between KDP(AHR) and A(ZPHI) demonstrated in section 3b and 2) the fact that the presented data were collected from one radar. Hence, the following analysis is based on internal polarimetry consistency.

The scatterplots A(CZPHI, C)–KDP(AHR) and A(CZPHI, AHR)–KDP(AHR) resulting from event E1 are shown in Fig. 6a. Observe that multiple A(CZPHI, C) estimates are smaller than those from A(CZPHI, AHR) as a consequence of selecting “small optimal” αΦDP(C) values. The correlation coefficient ρAK from A(CZPHI, C)–KDP(AHR) is equal to 0.78, while from A(CZPHI, AHR)–KDP(AHR) it is 0.98. Their corresponding standard deviations σAK with respect to A = αKDP are 0.28 and 0.05 dB km−1, respectively, where α values are given by αΦDP(AHR). In Fig. 6b, attenuation-corrected Z(CZPHI, C) and Z(CZPHI, AHR) are compared against Z(DP, AHR), where Z(DP, AHR) is obtained from KDP(AHR) and αΦDP(AHR). Their root-mean-square errors (RMSE) are equal to 1.67 and 0.10 dB, respectively, for Z(DP, AHR) ≥ 35 dBZ. This means that the attenuation-correctioned CZPHI method can lead to lower performance than the ZPHI method, comparing Fig. 6b with Fig. 4b. In this analysis, the RMSE was used instead of the mean bias to take into account the standard deviation of Z(CZPHI) estimates associated with the variability of α. The quantified errors used to evaluate the CZPHI method are summarized in Table 3.

Fig. 6.

(a) The A(CZPHI, C)–KDP(AHR) and A(CZPHI, AHR)–KDP(AHR) scatterplots resulting from event E1 at 1216 UTC are represented by the red and green dots, respectively. (b) As in (a), but for Z(CZPHI, C)–Z(DP, AHR) and Z(CZPHI, AHR)–Z(DP, AHR) scatterplots. In addition, the relation Z(CZPHI) = Z(DP) is indicated by the black line.

Fig. 6.

(a) The A(CZPHI, C)–KDP(AHR) and A(CZPHI, AHR)–KDP(AHR) scatterplots resulting from event E1 at 1216 UTC are represented by the red and green dots, respectively. (b) As in (a), but for Z(CZPHI, C)–Z(DP, AHR) and Z(CZPHI, AHR)–Z(DP, AHR) scatterplots. In addition, the relation Z(CZPHI) = Z(DP) is indicated by the black line.

Table 3.

Comparison results between A(CZPHI, C) and A(CZPHI, AHR) using KDP(AHR) as a reference for four storm events.

Comparison results between A(CZPHI, C) and A(CZPHI, AHR) using KDP(AHR) as a reference for four storm events.
Comparison results between A(CZPHI, C) and A(CZPHI, AHR) using KDP(AHR) as a reference for four storm events.

A similar analysis of A(CZPHI) is performed using KDP(C) as a reference instead of KDP(AHR) and the results are summarized next. The correlation coefficient between A(CZPHI, C) and KDP(C) is equal to 0.59 and smaller than those shown in Fig. 6a. This is because of the limited accuracy associated with KDP(C). The resulting RMSE between Z(CZPHI, C) and Z(DP, C) is equal to 0.82 and smaller than the case when Z(DP, AHR) is used as a reference. This is because Z(CZPHI, C) and Z(DP, C) are obtained from the same αΦDP(C) values, leading to similar attenuation correction results. Nonetheless, even if Z(DP, C) is set as a reference, their resulting RMSE is still larger than the one from Z(CZPHI, AHR)–Z(DP, AHR).

Attenuated z and zdr and attenuation-corrected Z(CZPHI, AHR) and ZDR fields from event E1 are displayed in Fig. 7. The Z(CZPHI, AHR) field restored attenuated z areas with PIA values up to 14 dB mostly on the north side of the storm cell, which is associated with rapid increments of ΨDP (see Fig. 3). A similar situation is observed by comparing the fields of zdr and ZDR, where enhanced areas of ZDR correspond to oblate raindrops. From the ZDR field, it seems that its lower bound is between −2 and −1 dB, which could be due to radar miscalibration rather than prolate-shaped particles, and therefore Z and ZDR fields may not represent calibrated measurements. Furthermore, the radial pattern presented in the zdr and ZDR fields may be associated with an azimuthal modulation as result of a metallic fence near the radar causing PBB effects (Giangrande and Ryzhkov 2005). Although such error sources may cause uncertainties on Z and ZDR, they do not seem to affect estimates of KDP and A by neither of the discussed methods and they do not influence the results of the presented analysis.

Fig. 7.

Event E1 at 1216 UTC. Fields of (a) z, (b) zdr, (c) Z(CZPHI,AHR), and (d) ZDR are illustrated. The black contours represent the 30-dBZ level.

Fig. 7.

Event E1 at 1216 UTC. Fields of (a) z, (b) zdr, (c) Z(CZPHI,AHR), and (d) ZDR are illustrated. The black contours represent the 30-dBZ level.

b. Event E2: Mini-supercell

The performance of the CZPHI method from event E2 at 1450 UTC is analyzed in a similar manner as for event E1 and the quantified errors are summarized in Table 3. The results show again that the CZPHI method performs better when α is given by αΦDP(AHR) instead of αΦDP(C). Nonetheless, event E2 shows specific signatures associated with the spatial distribution of raindrop size that can be used to illustrate the ability of selecting proper α values using the outcome of both KDP approaches.

The resulting Z(CZPHI, AHR) and ZDR fields at 1450 UTC, associated with PIA (PIADP) values up to 10 dB (1.6 dB), are shown in Fig. 8. In the Z(CZPHI, AHR) field, a significant gradient can be seen along the inflow edge of the storm (arrow 1), as well as a narrow echo appendage (arrow 2). An echo appendage typically curves in the presence of a mesocyclone process; however, this feature was not seen during the considered period. The ZDR field shows an area of significantly enhanced values along the inflow edge (arrow 3). This feature, commonly seen in supercell storms, is referred to as the ZDR arc signature as a result of possible size sorting processes (Kumjian and Ryzhkov 2008). The fields of KDP(C) and KDP(AHR) are also illustrated in Fig. 8. It is seen that the KDP(AHR) field retains the spatial variability of the storm better than the KDP(C) field while reducing negative KDP estimates. Note that both KDP fields show enhanced values along the inflow edge of the storm with values as high as 12° km−1 collocated with the ZDR arc. Estimates of KDP over the echo appendage, in both KDP fields, are not possible because of its narrow width of less than 1 km.

Fig. 8.

Event E2 at 1450 UTC. Fields of (a) Z(CZPHI, AHR), (b) ZDR, (c) KDP(C), and (d) KDP(AHR) are shown. The black contours indicate the 30-dBZ level, and the magenta contours in (c) show the −1° km−1 level. In addition, the Z gradient along the inflow edge, the Z narrow appendage, and the ZDR arc signatures are indicated by arrows 1–3, respectively. The low-level inflow in (a) is represented by the three arrows.

Fig. 8.

Event E2 at 1450 UTC. Fields of (a) Z(CZPHI, AHR), (b) ZDR, (c) KDP(C), and (d) KDP(AHR) are shown. The black contours indicate the 30-dBZ level, and the magenta contours in (c) show the −1° km−1 level. In addition, the Z gradient along the inflow edge, the Z narrow appendage, and the ZDR arc signatures are indicated by arrows 1–3, respectively. The low-level inflow in (a) is represented by the three arrows.

The selected values for αΦDP(C) and αΦDP(AHR) are given in Fig. 9 as a function of azimuth. Observe that the optimization of α using ΦDP(C) was possible only in three azimuthal profiles of the mini-supercell. This is because in multiple azimuthal profiles, the percentage of gates per profile with KDP(C) > 0° km−1 is less than 50%, which led to the selection of the constant α, avoiding suboptimal α values. This means that in those profiles, A is given by the ZPHI method, leading to a reasonable correlation ρAK as shown in Table 3. On the other hand, the optimization of α using ΦDP(AHR) occurred in multiple azimuthal profiles, resulting in values mostly larger than 0.34 dB (°)−1 in contrast to those resulting from ΦDP(C). According to Ryzhkov and Zrnić (1995) and Carey et al. (2000), such large values are expected in areas of big raindrops, which is consistent with the ZDR arc signature.

Fig. 9.

Event E2 at 1450 UTC. Selected values for α using ΦDP(C) and ΦDP(AHR) are given by the stemplots in red and green, respectively, as a function of azimuth.

Fig. 9.

Event E2 at 1450 UTC. Selected values for α using ΦDP(C) and ΦDP(AHR) are given by the stemplots in red and green, respectively, as a function of azimuth.

c. Event E3: Tornadic cell

This event was associated with a bow apex feature along the leading edge of the storm. According to Funk et al. (1999), cyclonic circulations can occur along or near the leading bow apex, which can produce tornadoes of F0–F3 intensity. For a detailed observation of event E3, only the southeast side of the Z(CZPHI, AHR), ZDR, KDP(C), and KDP(AHR) fields at 1955 UTC are shown in Fig. 10. The Z field shows a strong gradient along the leading edge (arrow 4), indicating a region of strong convergence and low-level inflow (white arrows). A bow apex attribute resulting possibly from a descending rear inflow jet (Weisman and Trapp 2003) is also noticeable (arrow 5). This feature seems to be associated with a rotation pattern in the form of a hook or weak-echo hole (Bluestein et al. 2007) (extended arrow 6) that caused wind and tornado damage as indicated in Table 1. It is also observed that the core of the weak-echo hole, whose inner diameter is approximately 0.75 km, is related to bounded weak ZDR and KDP values, located in the center of the white circles. It can be observed that KDP(AHR) preserves the storm structure better than KDP(C) because the AHR approach avoids a segmented KDP texture and negative KDP values, which are observed in the KDP(C) field. Maximum values of PIA and PIADP reached 18 and 3 dB, respectively, while fully attenuated areas (south side) occurred behind strong rain echoes associated with KDP values on the order of 10° km−1.

Fig. 10.

Event E3 at 1955 UTC. Fields of (a) Z(CZPHI, AHR), (b) ZDR, (c) KDP(C), and (d) KDP(AHR) are indicated. Various levels are shown: 40 dBZ (black contours), 1° km−1 (red contours), and −1° km−1 (magenta contours). The Z gradient along the inflow edge (arrow 4), the bow apex (arrow 5), and the weak-echo hole (arrow 6) are given in (a). The low-level inflow (white arrows) and the rotation pattern associated with the echo-weak hole (white circles) are also shown.

Fig. 10.

Event E3 at 1955 UTC. Fields of (a) Z(CZPHI, AHR), (b) ZDR, (c) KDP(C), and (d) KDP(AHR) are indicated. Various levels are shown: 40 dBZ (black contours), 1° km−1 (red contours), and −1° km−1 (magenta contours). The Z gradient along the inflow edge (arrow 4), the bow apex (arrow 5), and the weak-echo hole (arrow 6) are given in (a). The low-level inflow (white arrows) and the rotation pattern associated with the echo-weak hole (white circles) are also shown.

The resulting values of αΦDP(C) and αΦDP(AHR) as a function of the azimuthal sector [0°; 360°], not shown here, indicate that for most azimuthal profiles, α values are associated with a minimum error E, except in the azimuthal sector of [40°; 120°], where estimates of A were determined by the ZPHI method. This sector was related to light and uniform rain profiles, where ΔΦDP values are smaller than 10°. Optimal values of αΦDP(AHR) are predominantly found between 0.34 and 0.50 dB (°)−1. The absence of α> 0.50 dB (°)−1, in contrast to event E2, may indicate the lack of big drops present at this time. Selected values of αΦDP(C) are frequently smaller than or equal to 0.34 dB (°)−1 but in a few profiles they are equal to 0.1 or 0.6 dB (°)−1, possibly as a result of an inadequate optimization process. The resulting e¯min and σemin, together with ρAK, σAK, and RMSE are given in Table 3, showing that ΦDP(AHR) profiles lead to more reliable values of α and better estimates of A and Z.

d. Event E4: Irregular-shaped cell

In contrast to events E1–E3, E4 is mainly related to light rain with a few spots of moderate rain and it is not associated with any known reflectivity signatures. In addition, multiple radial paths with reflectivity echoes larger than 30 dBZ are mostly smaller than 5 km, in which PIA reached values of 2.5 dB, and only in few profiles it increased to 14 dB. The fields of Z(CZPHI, AHR), ZDR, KDP(C), and KDP(AHR) at 0558 UTC are shown in Fig. 11. Comparing the fields of Z and KDP, the KDP(AHR) field maintains the spatial structure of the storm better than KDP(C). It can be seen that the magnitudes of KDP(C) and KDP(AHR) are frequently smaller than 4° km−1, implying a slow incremental behavior of estimated ΦDP profiles. As such, only the azimuthal sectors [75°; 150°] (east side) and [250°; 280°] (west side) were associated with ΔΦDP> 10°. In both sectors, the optimization process was characterized by an inadequate performance because, in multiple azimuthal profiles, repetitive values equal to 0.1 dB (°)−1 were selected and the associated errors were larger than those found in E1–E3. In the remaining profiles, values of αΦDP(C) were smaller than 0.34 dB (°)−1, while values of αΦDP(AHR) were comparable to 0.34 dB (°)−1, indicating the absence of raindrops of considerable size. The results associated with the selection of α using ΦDP(C) and ΦDP(AHR) are indicated in Table 3, showing a decreased performance of the CZPHI method compared to the results of E1–E3.

Fig. 11.

As in Fig. 8, but for E4 at 0558 UTC.

Fig. 11.

As in Fig. 8, but for E4 at 0558 UTC.

5. Evaluation of δhv estimates

For each storm event, the preprocessed ΨDP (section 3a), the obtained KDP(AHR) fields, and A(CZPHI, AHR) fields were set as inputs to the δhv algorithm for its evaluation. As part of the δhv approach (step 1), a low-pass FIR filter specified by a 32-filter order and 1-km cutoff range scale was applied to the ΨDP field.

The estimated δhv fields resulting from storm events E1–E4 at 1216, 1450, 1955, and 0558 UTC, respectively, are shown in Fig. 12. In all events, it can be seen that the areas of δhv that are given by a uniform value correspond to the areas of Z smaller than the 30-dBZ level, which defines the shape of the described storm cells. Moreover, a spatial correlation between the δhv fields and their corresponding ZDR fields is observed, which confirms the correlation nature between δhv and ZDR (e.g., compare Figs. 12a and 7d). Such a spatial correlation is not exclusive to δhv and ZDR because a similar correlation is also observed between the fields of δhv, Z, and KDP, exemplifying the self-consistency relation (Scarchilli et al. 1996) between ZDR, Z, and KDP in a comparable manner.

Fig. 12.

The resulting fields of δhv from E1 to E4 at (a) 1216, (b) 1450, (c) 1955, and (d) 0558 UTC. The white contours indicate the 30-dBZ level, whereas the black contours in (c) represent the 40-dBZ level.

Fig. 12.

The resulting fields of δhv from E1 to E4 at (a) 1216, (b) 1450, (c) 1955, and (d) 0558 UTC. The white contours indicate the 30-dBZ level, whereas the black contours in (c) represent the 40-dBZ level.

The ability of the algorithm to capture the spatial variability of δhv is substantial. For example, in Fig. 12a, significant δhv values are more visible on the north side than on the south side of the storm cell, indicating the presence of Mie scattering. Another example of the spatial variability and consistency aspects of δhv