One of the most beneficial polarimetric variables may be the specific differential phase KDP because of its independence from power attenuation and radar miscalibration. However, conventional KDP estimation requires a substantial amount of range smoothing as a result of the noisy characteristic of the measured differential phase ΨDP. In addition, the backscatter differential phase δhv component of ΨDP, significant at C- and X-band frequency, may lead to inaccurate KDP estimates. In this work, an adaptive approach is proposed to obtain accurate KDP estimates in rain from noisy ΨDP, whose δhv is of significance, at range resolution scales. This approach uses existing relations between polarimetric variables in rain to filter δhv from ΨDP while maintaining its spatial variability. In addition, the standard deviation of the proposed KDP estimator is mathematically formulated for quality control. The adaptive approach is assessed using four storm events, associated with light and heavy rain, observed by a polarimetric X-band weather radar in the Netherlands. It is shown that this approach is able to retain the spatial variability of the storms at scales of the range resolution. Moreover, the performance of the proposed approach is compared with two different methods. The results of this comparison show that the proposed approach outperforms the other two methods in terms of the correlation between KDP and reflectivity, and KDP standard deviation reduction.
Weather observations are conventionally performed by single-polarization S- or C-band weather radars. Although these radars have substantially improved weather monitoring, researchers have found several limitations. For example, the spatial and temporal resolutions obtained from these conventional radars seem to be undesirable for the early detection of small but threatening features of convective weather as well as the detection of localized and heavy rainfall storms (Heinselman and Torres 2011; Schellart et al. 2012; Berne and Krajewski 2013). In contrast, single-polarization X-band weather radars are suited to obtain localized weather observations at resolutions higher than those of conventional radars. For example, a network of X-band weather radars in Hamburg, Germany, is used to observe precipitation at high resolutions (Lengfeld et al. 2016). Nonetheless, power attenuation and radar miscalibration may reduce the accuracy of single-polarization radar observables (Gourley et al. 2009).
One technique to potentially mitigate these issues is polarimetric radar technology (Doviak et al. 2000; Bringi and Chandrasekar 2001). For instance, in western Europe, polarimetric X-band weather radars are used to obtain high-resolution rainfall rates in order to cope with urban flooding (Ochoa-Rodriguez et al. 2015). In the United States, a network of X-band radars with polarimetric capabilities is used to improve the coverage of convective weather at low levels (McLaughlin et al. 2009), while in France a similar radar network is used to fill the gaps from operational radars in mountainous regions (Beck and Bousquet 2013). Moreover, in Germany, polarimetric Doppler X-band radars together with Doppler lidars are installed in international airports to provide local measurements of precipitation type and wind shear conditions, which are difficult to obtain from a national weather radar network (Ernsdorf et al. 2015).
In this work the polarimetric radar variable of interest is the specific differential phase KDP because it is independent of attenuation and miscalibration; therefore, it can improve rainfall-rate estimation compared to power-based measurements, usually in heavy rain (Wang and Chandrasekar 2010). However, the accuracy of KDP is not yet sufficient because the measured differential phase ΨDP could be significantly noisy in light and moderate rain. In addition, in moderate and heavy rain, ΨDP can include a nonnegligible scattering component resulting from the Mie scattering region, which is known as the backscatter differential phase δhv (Matrosov et al. 2002; Trömel et al. 2013). Thus, accurate estimation of KDP is necessary in order to unleash the potential of polarimetric weather radars.
Literature review reveals a large and continuous study to estimate KDP and for simplicity it is divided into two groups. For the first group and in situations where δhv can be neglected (e.g., at S-band frequency or light rain), straightforward approaches based on autoregressive average models are applied to smooth ΨDP (Bringi and Chandrasekar 2001; Matrosov et al. 2006; Vulpiani et al. 2012). For the same group but in cases where δhv is of significance, Hubbert and Bringi (1995) introduced an iterative filtering approach to smooth ΨDP and filter δhv. A common problem in this group is that KDP is estimated with inadequate spatial resolutions that could result in an underestimation of KDP peaks and therefore lead to an inaccurate phase-based rainfall estimation (Ryzhkov and Zrnić 1996). This limitation was reduced by Wang and Chandrasekar (2009), who developed an algorithm to filter δhv and to control the smoothing degree on ΨDP while maintaining its spatial resolution. For the second group, the KDP approaches included polarimetric relations in rain, such as the self-consistency (SC) relation, which formulates a dependency between KDP, reflectivity Z, and differential reflectivity ZDR (Scarchilli et al. 1996; Goddard et al. 1994; Gorgucci et al. 1992), and the relation between δhv and ZDR, δhv–ZDR (Scarchilli et al. 1993; Testud et al. 2000). Otto and Russchenberg (2011), were able to estimate KDP at range resolution scales using both relations, while Schneebeli and Berne (2012) included the δhv–ZDR relation in their Kalman filter approach. Recently, Giangrande et al. (2013) introduced a linear programming method that includes Z measurements, whereas Huang (2015) used the SC relation to estimate KDP at S-, C-, and X-band frequencies. A disadvantage of using polarimetric relations is that uncertainties on Z and ZDR measurements might reduce the performance of these approaches. Last, approaches in both groups could be associated with significant errors when KDP is estimated at range resolution scales (Grazioli et al. 2014; Hu et al. 2015).
In this work an adaptive approach that includes polarimetric relations is presented to estimate accurate KDP from ΨDP in rain, whose δhv is of significance, at high spatial resolution. This paper is a follow-up of Otto and Russchenberg (2011) and is organized as follows. Two KDP methods—given by Hubbert and Bringi (1995), first group, and Otto and Russchenberg (2011), second group—are shortly described in section 2. They will be used for comparison with the proposed technique because (i) the method of Hubbert and Bringi (1995) is widely accepted for KDP estimation at C- and X-band frequency and (ii) the focus of this work is to improve the technique introduced by Otto and Russchenberg (2011). In section 3, the adaptive high-resolution approach is introduced to estimate KDP and model its standard deviation. To demonstrate the capability of this approach in terms of δhv contamination and spatial resolution, one storm event is analyzed in section 4. In section 5, the performance of the adaptive approach is compared with those from section 2, using four storm events. In section 6, conclusions are drawn. Finally, the standard deviation of the proposed KDP estimator is derived in appendix A, while the filter required by Hubbert and Bringi (1995) is designed in appendix B.
2. Specific differential phase: Background and estimation
In polarimetric weather radars, the difference between the horizontal and vertical polarization phases is defined as ΨDP (°). A conceptual model for a ΨDP profile is expressed as
where ΦDP(r) (°) represents the cumulative propagation phase shift along its course, while δhv(r) (°) indicates local backscattering phase shifts manifested as “blips” or “bumps.” Random noise is represented by ε (°) and the range by r (km). Depending on the weather environment, the standard deviation error of ΨDP, σP (°), varies between 2° and 5° (Bringi and Chandrasekar 2001).
The one-way KDP (° km−1) is half the derivative of ΦDP:
however, the estimation of accurate KDP is not straightforward. For rainfall-rate applications, KDP should be estimated such that the normalized standard error (NSE) of rainfall rate is less than 20% (Gorgucci et al. 1999; Bringi and Chandrasekar 2001). For instance, assuming a standard deviation of ΦDP is equal to 3° and estimating KDP as in Eq. (2) over a pathlength of 2 km, the standard deviation of KDP, σK, would be 1° km−1. If the rainfall rate and KDP are given by a power-law relation with a power coefficient of 0.8 (Ryzhkov et al. 2005; Wang and Chandrasekar 2010), then the value of σK is theoretically sufficient for KDP larger than 4° km−1. However, for KDP values smaller than 4° km−1, σK values are required to decrease accordingly, which can be achieved by increasing the pathlength with the sacrifice of spatial resolution.
a. Conventional method
Hubbert and Bringi (1995) introduced an iterative range filtering technique based on two steps. In the first step a low-pass filter is designed such that fluctuations resulting from δhv and ε at scales of the radar resolution [Δr (km)] are eliminated from the ΨDP profile. The magnitude response H of this filter is given in the range domain by defining its Nyquist frequency as 1/(2Δr). In this domain, H is specified by the filter order and Δr/rc, where rc (km) is the required cutoff scale such that rc > Δr. This means that the filter will maintain (or “pass”) spatial variations of ΨDP(r) at scales larger than rc. However, spatial variations smaller than rc will not be effectively suppressed because H is characterized by a transition bandwidth from the “pass” band to the “stop” band. This transition can be faster if the order of the filter is heavily increased. However, a high-order filter will strongly smooth the spatial variability of ΨDP, leading to a coarse spatial resolution of KDP. Thus, the order of the filter is selected such that H will suppress spatial variations at scales slightly larger than Δr without compromising the spatial resolution of KDP. This filter is referred to as a “light” filter that will lead to reduced σP and thereby σK.
In the second step the light filter is applied several times to eliminate extended δhv fluctuations, up to rc scales, as it would result from a “heavy” filter but without excessive smoothing. This process begins by filtering ΨDP, resulting in a first estimation of ΦDP . The absolute difference between ΨDP and profiles at each range gate are employed to generate a modified ΨDP profile (): if the absolute difference is larger than a threshold τ, then = holds; otherwise, = ΨDP(r) holds. This threshold is predefined from the interval [1.25; 2]σP. The process is iterated until from two consecutive iterations show insignificant changes. Next, from the last iteration is filtered one more time to obtain ΦDP and thereby KDP according to Eq. (2).
Even though this is an elegant approach to estimate KDP in real time, the following issues limit its purpose. First, spatial fluctuations larger than rc will not be completely eliminated by the iterative step. Second, its performance is sensitive to the value of τ. For example, τ = 2σP will lead to a more correlated to ΨDP than to , which might not be sufficient to eliminate unwanted fluctuations. Third, if rc is increased or Δr is decreased, then the order of the filter should be increased (Wang and Chandrasekar 2009), which may exacerbate these issues. In summary, in the design of the filter and the selection of the threshold, there is a trade-off between the smoothing extent and the spatial resolution required to estimate KDP with small standard deviation and reduced bias.
b. High-resolution method
In contrast to the conventional approach, Otto and Russchenberg (2011) included δhv–ZDR and SC relations in rain to estimate KDP at X-band frequency. The first relation is represented by δhv = and the second by
where Z and ZDR are given in dBZ and decibels (dB), respectively, while coefficients e1, e2, c1, c2, and c3 establish the average fit from a set of drop size distributions (DSDs), drop shape models, and temperatures. The specific differential phase in Eq. (3) is indicated by to distinguish between KDP from the SC relation and KDP from the high-resolution approach. Otto and Russchenberg (2011) used a normalized gamma distribution to model DSDs. To represent rain variability, 1500 DSDs were modeled by varying the parameters of the distribution (i.e., median volume diameter, concentration, and shape parameter). In addition, a hybrid drop shape model that consists of three models and temperatures in the range of 1°–25°C was considered to express a wide range of possibilities. The coefficients for the δhv–ZDR relation are e1 = 1 and e2 = 1.8, while for the SC relation they are c1 = 1.37 × 10−3, c2 = 0.68, and c3 = −0.042.
The measurements for Z and ZDR are corrected for attenuation and differential attenuation, respectively, according to Bringi et al. (1990). The path-integrated attenuation (PIA) in reflectivity and in differential reflectivity (PIADP) is given by 0.34ΔΨDP and 0.05ΔΨDP, respectively, where both coefficients are obtained from scattering simulation. For a ΨDP profile represented by Eq. (1), the difference of ΨDP in a path can be approximated as ΔΨDP = ΔΦDP + Δδhv, where ΔΨDP = ΨDP(b) − ΨDP(a) and b > a. The length of could be as small as Δr or as large as the maximum unambiguous range. To identify whether Δδhv is negligible, the δhv–ZDR relation is used in the following assumption: if |ZDR(b) − ZDR(a)| is smaller than 0.3 dB, then δhv(b) − δhv(a) ≈ 0°. If that case is satisfied, then ΔΨDP is retained for further processes; otherwise, ΔΨDP is discarded. Multiple ΔΨDP samples associated with negligible Δδhv are obtained by considering more paths. These samples are weighted by Z and ZDR using the SC relation to obtain ΔΨDP samples at Δr scale. For instance, the weight w at gate i within path is
Weighted ΔΨDP samples are used according to Eq. (2) to obtain multiple KDP samples. Finally, KDP(i) is estimated by the arithmetic mean of these samples.
Otto and Russchenberg (2011) demonstrated their method using one rainfall event, showing a visual consistency between KDP and Z as well as between δhv and ZDR. Nevertheless, fewer consistent results were observed for KDP values smaller than 4° km−1. In addition, estimates of KDP were associated with values of σK as high as 3° km−1. Moreover, the sensitivity of the SC relation to DSD, drop shape, and temperature variability and its impact on estimated weights were not discussed. In addition, the effect of uncertain measurements of Z and ZDR caused by, for example, attenuation, partial beam blockage (PBB), and miscalibration on the performance of this approach remains an open question. Besides its limited validation, another matter is the computational time required to weigh a significant amount of ΔΨDP samples at each range gate, which might be inefficient for operational purposes.
3. Adaptive high-resolution approach
The presented approach is an improved version of the high-resolution method in order to address limitations associated with the conventional and high-resolution methods, mainly the low accuracy of KDP in light rain, contamination resulting from unfiltered δhv, sensitivity to Z and ZDR measurements, and computational time. The inputs are radial measurements of ΨDP, Z, and ZDR in rain. In addition, a predetermined length interval [Lmin; Lmax] is required to control the selection of pathlengths. This interval is assumed to be defined by a user; however, possible values are discussed in section 3e. The adaptive approach consists of three processes: preprocessing, pathlength selection, and KDP estimation. A flowchart to estimate KDP, among other variables, for a given radial profile, is presented in Fig. 1.
To correct Z and ZDR profiles for attenuation and differential attenuation, respectively, Otto and Russchenberg (2011) used noisy ΔΨDP instead of ΔΦDP, which may decrease the accuracy of the method given by Bringi et al. (1990). In contrast, in the adaptive high-resolution approach, a linear regression fit over a 3-km window is applied to a ΨDP profile, resulting in a ΦDP profile (). Thus, attenuation-corrected reflectivity (Zt) and differential reflectivity () are given as and , respectively, where z (dBZ) and zdr (dB) represent attenuated and differential attenuated measurements, respectively. This attenuation correction method might be sensitive to measurements errors, constant coefficients, and δhv contamination. However, Gorgucci and Chandrasekar (2005) studied the method of Bringi et al. (1990) at X-band frequency and showed that this method performs fairly well with only a slight degradation of the average error for attenuation correction.
In this approach an estimate of the standard deviation of profile, σZDR (dB), is required by the pathlength selection process. A moving window of five gates is applied to the profile so that local σZDR samples are obtained. Then, σZDR is estimated by averaging these σZDR samples. The estimation of KDP is achieved gate by gate, starting from ranges near the radar and continuing downrange. Assuming that the first gate with measurements of rain is located at ri, the estimation of KDP begins at gate i.
b. Pathlength selection
In the high-resolution technique, ΔΨDP samples were obtained from ΨDP using paths of any possible length. However, KDP results were associated with high values of σK and significant computational time. In this work a pathlength L (km) for gate i is selected from [Lmin; Lmax] such that a theoretical σK is minimized. The formulation of a theoretical σK is shown in section 3e but now let the theoretical σK be a function of parameters L and M, where M represents the number of ΔΨDP samples with negligible Δδhv. To identify negligible Δδhv, the high-resolution technique used a fixed threshold to constrain . However, a fixed threshold might not capture the possible variability of ZDR within the storm. In this work the condition to identify negligible Δδhv is given by
Equation (5) can be considered to be independent of any parameterized δhv–ZDR relation because this relation is not used quantitatively. Instead, Eq. (5) relies on the existing correlation between δhv and ZDR. Furthermore, the sudden variability in microphysics is taken into consideration by using σZDR rather than an arbitrary threshold. Issues such as ZDR miscalibration are mitigated by the estimation of . Equation (5) is referred to as the ΔΨDP filter condition.
The pathlength selection starts with L = Lmin. For simplicity, a pathlength is of the form L = nΔr, where n is an integer larger than 1. Then, a range interval centered at gate i is defined as [ri − L; ri + L]. This range interval is used to limit the extent of and to obtain multiple samples. These samples are achieved by shifting a path of length L, within the interval starting at ri − L, n times with steps of Δr. In this manner, (n + 1)− samples are obtained. Next, M is calculated by counting the number of samples that satisfy Eq. (5). Note that M ≤ (n + 1). From L and M, a σK value is determined. To have a set of σK values, the same procedure is repeated for the next value of L until L = Lmax. Finally, the pathlength that leads to the minimum σK is selected and is represented by .
Repeating a similar procedure but with [i.e., shifting a path of length and using Eq. (5)], M–ΔΨDP samples with negligible Δδhv are retrieved from ΨDP to estimate KDP(i). The remaining ΔΨDP samples are discarded to avoid bias on KDP(i). For the next steps, only ΨDP, Zt, and profiles in the interval are used.
To estimate KDP(i), M–ΔΨDP samples should be downscaled from to Δr scales. A downscaling weight w(i) was suggested by the high-resolution method according to Eq. (4). In contrast, in the adaptive high-resolution approach, a different formulation of w(i) is proposed in order to reduce its sensitivity to possible sources of uncertainty that were discussed in section 2b, mainly the sensitivity of the SC relation to rain variability, radar miscalibration, and PBB effects. Moreover, this formulation allows us to study statistics of w(i) and KDP(i) for quality control purposes.
Consider a theoretical ΔΨDP ≥ 0° from a path of length L. For gate i in the interval [a + 1; b], the downscaling weight w(i) is expressed as a factor that weighs ΔΨDP such that ΔΨDP(i) = w(i)ΔΨDP, where ΔΨDP(i) represents the difference ΨDP(i) − ΨDP(i − 1) (i.e., at Δr scale). For derivation purposes, let w(i) be bounded by the interval [0; 1] and 1. Using Eqs. (1) and (2), ΔΨDP(i) and ΔΨDP are expressed in terms of KDP and w(i) as
Both KDP(i) and KDP are estimated using the SC relation according to Eq. (3) at scales of Δr and L, respectively. In the numerator of Eq. (6), adjacent δhv values are assumed to be similar, so δhv(i) − δhv(i − 1) is approximately 0°. In the denominator, assuming that ΔΨDP satisfies Eq. (5), the difference δhv(b) − δhv(a) is negligible. In this manner, w(i) is formulated as
where the SC relation is used two times in contrast to the high-resolution method, which is used (n + 2) times [see Eq. (4) and replace b with a + n]. In this way, added errors associated with the SC relation are reduced.
To downscale M–ΔΨDP samples (i.e., with j = 1, 2, … M) associated with M– paths of length , Eq. (7) is used and the jth weight is given as
where and represent the arithmetic mean of and values in path , respectively. Repeating Eq. (8) over the remaining paths, M– samples are obtained and KDP(i) is estimated as
Once KDP(i) is estimated, the pathlength selection and KDP estimation processes are applied to gate i = i + 1 until the last gate measured in rain. Hence, a KDP profile is obtained as well as associated and M profiles.
The KDP estimator is a function of variables and , which result from the ΔΨDP filter condition and ΔΨDP downscaling, respectively. Therefore, it is important to discuss errors associated with both variables. For this purpose, is expressed as , where εδ indicates possible errors from neglecting Δδhv. Using the scattering simulation introduced in section 2b and setting σZDR equal to 0.2 dB, the estimated mean and standard deviation of εδ are 0.04° and 0.6°, respectively. The uncertainty of depends on the SC relation in rain. Trabal et al. (2014) demonstrated that the coefficients of the SC relation shown in Eq. (3) are sensitive to temperature variability, while DSD and drop shape variabilities are well represented by a normalized gamma distribution and a hybrid drop shape model. Similar findings were reported by Gourley et al. (2009) and Adachi et al. (2015). Although is independent of c1, any possible sensitivity to c2 and c3 is modulated by the difference and , respectively (i.e., it depends on the spatial variability of Zt and within path instead of their absolute values). For example, in a uniform path, might be constant and equal to . Moreover, is independent of constant biases in Zt and within as well as radar miscalibration. This independence could reduce the impact of biases on Zt and areas caused by PBB. For simplicity the estimated weight is rewritten as = , where is referred to as the SC ratio [see Eqs. (7) and (8)].
The uncertainty of can be quantified by its NSE. Scarchilli et al. (1996) derived an expression for NSE of , hereinafter NSE(K), that is a function of c2, c3, and variances of Zt and . Using this expression and basic properties of the variance, the NSE of is given by NSE = NSE(K). For example, setting n equal to 10 (i.e., L is 10 times Δr) and using values of c2 and c3 given in section 2b and conventional accuracies of 1 and 0.2 dB for Zt and , respectively, NSE results in 15.7%. This analysis can be used as guidance to identify which elements associated with the ΔΨDP filter condition and ΔΨDP downscaling can lead to an incorrect estimation of KDP.
The uncertainty of is measured by its standard deviation as
Equation (10) is referred to as the actual σK. In addition, the NSE of , hereinafter , is given by 100%. Both actual σK and profiles are added to the output of Fig. 1, which can be used for quality control purposes. In a similar manner but for the M– samples, profiles of their actual mean (μα), standard deviation (σα), and NSE () are also obtained.
The uncertainty of can be controlled by modeling its actual σK. Therefore, a theoretical σK is derived from Eq. (9), where is replaced by . Values of Δr, L, and M are assumed to be given, while and are defined as random variables. Using properties of the variance that involve the sum and product of random variables, a theoretical σK is approximated as (for a detailed derivation, refer to appendix A)
where and are the variance of ΨDP and εδ, respectively. The value of μα depends on the spatial variability of and . For example, μα may be near 1 in stratiform rain, but it can be between 2 and 5 in convective rain. The value of σP also depends on rain type, while σε is assumed equal to 0.6° as given in section 3d. To illustrate theoretical σK, let us assume that μα and σP are equal to 3 and 3°, respectively. A similar value for σP was given by Lim et al. (2013).
Theoretical σK curves as a function of M and for different combinations of L are presented in Fig. 2. Although σK in Eq. (11) is independent of Δr, the maximum value of M for a fixed pathlength L is given by (n + 1), which is equivalent to L/Δr + 1. For Δr equal to 0.03 km, Fig. 2a shows that if L is near 1 km, then it is expected to obtain σK values larger than 1° km−1. However, σK values smaller than 0.5° km−1 are expected if L is equal to 2 km and M is larger than 40. In terms of probability, M larger than 40 indicates that at least 60% of the total number of paths (n + 1) satisfies Eq. (5). For L larger than 3 km, σK curves continue for values of M beyond 100 but are not shown here. Note that the gap between consecutive σK curves decreases as L increases (e.g., σK curves corresponding to 5 and 6 km are almost identical). Therefore, it is recommended to avoid large values of L (i.e., lengthy paths) associated with a small reduction of σK.
Recall that for a given value of L, M is determined only after evaluating Eq. (5). To search for the combination of L and M that leads to the smallest theoretical σK, a pathlength interval [Lmin; Lmax] is considered instead of a single pathlength. For example, a length interval equal to [3;5] km could be defined from Fig. 2a as an input to the adaptive high-resolution approach. In that case, σK is expected to be small with a sufficient number of paths M. Similar σK curves can be produced for a coarser spatial resolution. For instance, if Δr = 0.25 km, then a set of L values ranging from 4 to 14 km would be used to generate σK curves as shown in Fig. 2b and an interval equal to [6;10] km would be defined. Note that for larger Δr, a smaller number of paths M might lead to reasonable σK values (cf. Figs. 2a and 2b).
In this σK modeling, it is assumed that a user predefines a theoretical value for μα and σP according to the storm type, and sets [Lmin; Lmax]. For example, in a uniform Z field (i.e., low μα and σP) or a more variable Z field (i.e., high μα and σP), [Lmin; Lmax] can be decreased or increased according to theoretical σK curves. In summary, setting [Lmin; Lmax] allows us to avoid high values of actual σK as well as unnecessary lengthy paths associated with an increased dependence on the SC ratio and a large computational time.
4. Analysis of the adaptive high-resolution approach
a. Data settings
The polarimetric X-band weather radar the International Research Centre for Telecommunications and Radar (IRCTR) drizzle radar (IDRA) is a frequency-modulated continuous wave system whose operational range is 15.3 km with a resolution of 0.03 km (Figueras i Ventura 2009). IDRA is located at the Cabauw Experimental Site for Atmospheric Research (CESAR) observatory in the Netherlands at a height of 213 m from ground level (Leijnse et al. 2010). It scans at a fixed elevation of 0.50° and rotates the antenna over 360° in 1 min with a beamwidth of 1.8°. Clutter echoes are removed by a filter based on spectral polarimetric processing (Unal 2009). Moreover, the measured linear depolarization ratio LDR is used to filter areas that include particles other than rain and have low SNR, such that range gates with LDR larger than −18 dB are removed (Bringi and Chandrasekar 2001). This simple filtering procedure should be extended in the case of an automatic algorithm. The unwrapping of differential phase profiles is performed by detecting a differential phase jump between two adjacent gates, equal to 80% of the maximum differential phase of 180°.
A convective storm event was observed by IDRA on 10 September 2011. The plan position indicators (PPIs) of z, zdr, and ΨDP are shown in Fig. 3. Attenuation-affected areas can be identified behind strong reflectivity echoes. The radial pattern observed in the fields of zdr and ΨDP is probably due to a metallic fence that surrounds the IDRA platform. This introduces an azimuthal-dependence bias in the zdr field in a similar way as PBB effects (Giangrande and Ryzhkov 2005). In addition, z might suffer from radar miscalibration as reported by Otto and Russchenberg (2011). Here, we use the opportunity to study the impact of rain attenuation and biases associated with PBB and miscalibration on the adaptive high-resolution approach.
The proposed KDP approach is analyzed using the storm event observed by IDRA. Besides the z, zdr, and ΨDP fields, a pathlength interval is required. For this requirement, Fig. 2a can be used for guidance because the theoretical σK curves were built for the same range resolution as for IDRA and with theoretical μα and σP equal to 3 and 3°, respectively, which express the spatial variability of the observed z and ΨDP fields. Thus, from Fig. 2a, a length interval of [3;5] km, which is associated with σK < 0.5° km−1 for M > 20, is selected. This interval is used by the pathlength selection process for all radials. For simplicity, is indicated by L(i) and is determined by minimizing instead of σK, as shown by Eq. (11), because theoretical μα, σP, and σε remain constant in all radials. The coefficients c2 and c3 given in section 2b are used by the KDP estimation process.
To study the performance of the KDP approach in terms of spatial resolution and δhv filtering, three variations on applying the ΔΨDP filter condition and ΔΨDP downscaling, given by Eqs. (5) and (9) respectively, are defined and indicated in Table 1. In case I, KDP is estimated using the ΔΨDP filter condition without downscaling (i.e., = Δr/L and the SC ratio is equal to 1). This case is denoted by Δδhv ≈ 0° and α ≈ 1. In case II, the ΔΨDP filter condition is not used but is downscaled. This case is expressed by Δδhv ≠ 0° and α ≠ 1. In case III, the ΔΨDP filter condition is applied and is downscaled, and this case is represented by Δδhv ≈ 0° and α ≠ 1. Note that case III follows the proposed KDP approach and that cases I and II are defined only for analysis purposes.
After a KDP profile is obtained, a ΦDP profile is reconstructed by integrating Eq. (2) and δhv is determined as δhv = ΨDP − ΦDP. Moreover, Z and ZDR are obtained by correcting z and zdr in a similar manner as in section 3a but replacing with ΔΦDP. This correction method could be improved by using more sophisticated techniques, such as those given by Park et al. (2005a) and Snyder et al. (2010). However, attenuation correction for z and zdr is beyond the scope of this work, as our goal is to assess the performance of the KDP approach.
c. and results
1) For a PPI radial
The KDP approach specified by cases I–III is applied to the azimuthal radial of 213° and its results are shown in Fig. 4. The downscaling aspect of KDP is examined by comparing the ΨDP, ΦDP (case I), and ΦDP (case III) profiles as shown in Fig. 4a. Observe that the total ΔΦDP for cases I and III is equal to 45°. However, ΦDP from case III captures the spatial variability and rapid increments of ΨDP better than ΦDP from case I. This can be seen by their corresponding KDP profiles, which are also shown in Fig. 4a but they are shifted by −10° km−1, where two KDP peaks (both of approximately 10° km−1) from case III correspond to fast increments of ΨDP located downrange in convective areas. Observe that the ΨDP profile includes a δhv “bump” in the range 2–4 km. To analyze the δhv contamination aspect of KDP, ΦDP from cases II and III are shown in Fig. 4b as well as their corresponding δhv profiles shifted by −10°. This δhv bump of 2-km length does not show an impact on ΦDP (case II) because its length is smaller than Lmin = 3 km. Also, ΨDP values outside the bump are similar and therefore most ΔΨDP samples have values near 0°. In summary, KDP estimation is not affected by this δhv bump in both cases, II and III. Note that in the range 7–11 km, ΨDP increases rapidly, which probably means that raindrops of moderate to large size are present and thus significant δhv values are expected. However, both ΦDP profiles are similar. This similarity may be due to δhv values hardly varying in this range and thus Δδhv samples do not impact ΔΨDP samples. Such a feature can be seen in δhv (case III), where it shows a slight spatial variability. Further, the estimation of δhv depends on ΦDP estimation, which can include accumulated KDP errors (e.g., toward the end of the range). A rigorous estimation and analysis of δhv are beyond the scope of this work; the focus here is on KDP estimation.
Figure 4c shows attenuated z, corrected Zt, and corrected Z profiles—the last profile being associated with case III. The correction of z is evident toward the convective range 7–11 km, where PIA reaches 15 dB at 11 km. Note that the Zt profile shows values slightly larger than those of the Z profile because of unfiltered δhv. In a similar manner, the correction of the zdr profile is shown in Fig. 4d, where PIADP equals 2.3 dB at 11 km. Note that the δhv and ZDR profiles show a correlated behavior as expected from the δhv–ZDR relation.
2) For a full PPI
The results from applying the KDP approach, specified by case III, to all radials of Fig. 3 are shown in Fig. 5. The field of Z is plotted in Fig. 5a and the field of L, selected from the interval [3; 5] km, in Fig. 5b. The spatial variability of L exhibits an adaptive performance with the purpose of minimizing σK. The KDP and ΦDP fields are shown in Figs. 5c and 5d, respectively. It can be seen that the KDP field maintains the structure and resolution of the storm illustrated by the Z field, whereas the ΦDP field displays the propagation phase component of the ΨDP field depicted in Fig. 3c. Note that the KDP field presents some gaps in areas of measured Z (e.g., at approximately 8 km south from IDRA). In these areas, M is equal to 0, which means that for any bounded value of L, none of the (n + 1)–ΔΨDP samples satisfies Eq. (5). Such an issue could be avoided if, for instance, KDP is estimated using case II instead of case III or interpolation algorithms are used. Figure 5e shows the actual σK, whose arithmetic mean is equal to 0.1° km−1. However, values as high as 1° km−1 are visible near convective edges. This increase in σK is partly due to a reduced number of ΔΨDP samples that satisfy Eq. (5). The actual μα are represented in Fig. 5f, whose values are mostly between 0 and 5. Thus, setting theoretical μα equal to 3 in this convective storm is a reasonable predefinition. The field of μα also shows an adaptive characteristic of the KDP approach as it handles the spatial variability of ΨDP. In a similar manner, KDP is estimated using cases I and II.
To study the impact of unfiltered δhv on the standard deviation of KDP, actual σK resulting from cases II and III over all radials is displayed in Fig. 6 as a function of its corresponding number of paths M, with M > 1. Because of the large number of σK and M samples, they are plotted as 2D histograms for better visualization. Note that the σK–M histogram from case II shows a very small dependence on M as opposed to Eq. (11). This holds for values of M up to 167 (i.e., the nearest integer less than Lmax/Δr + 1). Such behavior occurs because the ΔΨDP filter condition is not applied in case II and thereafter adjacent paths are employed to obtain ΔΨDP samples, which leads to an increased correlation coefficient between these samples. In a hypothetical situation with a correlation coefficient equal to 1, theoretical σK is no longer a function of M [see Eq. (A2) in appendix A]. In addition to these adjacent samples, unfiltered Δδhv, and thereby δhv, compromises the estimation of KDP and increases the variability of actual σK. In contrast, the 2D histogram from case III shows a dependence on M because σK decreases when M increases as expressed by the theoretical σK in Eq. (11). Two theoretical σK curves, for the same range resolution of IDRA, are also plotted in Fig. 6 to compare theoretical σK with actual σK from case III. The upper curve is set with L = 3 km (i.e., Lmin) and μα = 5 and the lower curve with L = 5 km (i.e., Lmax) and μα = 0.5, assuming the same theoretical values for σP and σε as given in section 3e. For a fair comparison between actual σK near 0° km−1 and the lower σK curve, only actual σK values slightly larger than 0° km−1—for example, ≥0.05° km−1—are considered. In this comparison, 91% of the σK–M scatters are in between both curves, while only 2% are located above the upper curve.
Another manner to study the downscaling and δhv contamination aspects of the KDP approach is by examining the consistency between Z and KDP (Park et al. 2005b). For this purpose, Fig. 7a compares the Z–KDP histograms from cases I–III. In case I KDP is estimated at coarse resolution and its values are smaller than 8° km−1, while in cases II and III KDP is estimated at Δr scales and KDP values can be as high as 12° km−1. However, in case II the 2D histogram shows multiple negative outliers because of unfiltered δhv, resulting in underestimated and overestimated KDP values. Among these three cases, case III provides the best consistency because of the application of the ΔΨDP filter condition and ΔΨDP downscaling specified by the proposed KDP approach.
For evaluation purposes, Fig. 7b shows the Z–KDP histogram from case III, the Z–KDP scatterplot using scattering simulation from section 2b, and its theoretical fit Z–KDP relation given by KDP = 8.7 × 10(0.69Z/10)−4. Note that simulated Z values were shifted by −8 dB in order to match those from case III, which could be due to incorrect attenuation correction and/or bias associated with PBB and miscalibration. In contrast, the KDP axis shows a noticeable agreement between simulation and estimation. As a first step to analyze the discrepancy in the Z axis, a similar histogram is shown but with attenuated z instead of corrected Z, keeping estimated KDP. From both plots it is clear that attenuation is not the major reason for this inconsistency but rather PBB and miscalibration.
d. in solid or mixed precipitation
As part of the presented analysis, KDP estimation at X-band frequency over areas of solid and melting hydrometeors, such as graupel, hail, and snow, are shortly discussed. Because the shape and orientation response of particles are strongly related to their dielectric response, polarimetric signatures of solid hydrometeors is reduced because their dielectric constant factor is 20% or less that of raindrops. For example, Dolan and Rutledge (2009) and Snyder et al. (2010) simulated KDP values at X band for solid and melting particles, showing a limited range of −1 to 1 (° km−1) except for melting graupel, which can be between −2 to 7 (° km−1). In addition, values of δhv from solid hydrometeors are small except for large nonspherical hail or melting hail, in which δhv can be in the order of 4°–7° (Trömel et al. 2013). Moreover, Schneebeli et al. (2014) calculated KDP in snow using a Kalman filter–based approach and found similar results when KDP is estimated by the conventional technique. In this context the spatial variability of KDP in nonrain regions may be less significant than in rain regions. Thus, the conventional approach or an autoregression-based model can be considered to complement KDP estimates in nonrain regions as demonstrated by Lim et al. (2013). Alternately, the adaptive high-resolution approach can be also used by setting equal to ΔrL−1 and M equal to n + 1, which is similar to case I but without the ΔΨDP filter condition, at expenses of low-resolution and possibly δhv infiltration. In this scenario the theoretical σK is simplified to assuming σP = 3°, μα = 1, and M = n. For instance, for Δr values of 0.03, 0.25, and 1 km, L can be set to 3, 6, and 9 km, respectively, in order to obtain a theoretical σK equal to 0.07° km−1. However, further research is required to test the suggested alternative.
5. Assessment of the adaptive high-resolution approach
In this section the performance of the proposed KDP approach specified by case III is compared with the conventional and high-resolution (HR) techniques. For this purpose, four storm events (E1–E4) observed by IDRA are described in Table 2. Although only observations of E3 at 1950 UTC were shown in section 4, the other events also display patterns related to attenuation, PBB, and miscalibration. In the conventional technique, the filter is designed using a 36th-order filter with rc = 1 km and τ = 1.5σP. More details on the filter design are included in appendix B.
a. During 1 min
The corresponding times for E1–E4 are 2151, 2225, 1950, and 0550 UTC, respectively. The Z–KDP histograms resulting from the three KDP approaches applied to each 1-min event (i.e., one PPI) are shown in Fig. 8. The Z–KDP scatters from the conventional technique are significantly spread because KDP is estimated at coarse resolution and δhv is not properly filtered, which leads to negative and positive KDP bias. In contrast, results from the HR method show more condensed relations. However, for Z values smaller than approximately 40 dBZ, multiple outliers are noticeable. Those outliers are substantially eliminated by the adaptive HR approach, which exhibits an enhanced consistency for weak and strong Z. To quantify the consistency of the results, the correlation coefficient between Z and KDP, ρZ,K, obtained from each approach is given by the second, third, and fourth columns of Table 3. From this quantification, the adaptive HR approach outperforms the other two techniques. For reference purposes, the ρZ,K from the simulated Z–KDP shown in Fig. 7b was also estimated and is equal to 0.75, which is similar to those resulted from the adaptive HR technique. Note that ρZ,K values resulting from simulation or observations can change according to the range of Z and KDP values because the theoretical Z–KDP relation is nonlinear and thereby ρZ,K may be used as a relative quantity.
The results of Z and KDP from the conventional and adaptive HR methods, presented in Fig. 8, show a similar discrepancy in the Z axis as indicated in Fig. 7b. Although the degree of discrepancy is not the same in all events, the Z values reached by the conventional technique are in the same order as those from the adaptive HR approach. This indicates that Z is most likely biased and that the estimation of KDP by the adaptive HR approach may not be affected by attenuation and biases associated with PBB and/or miscalibration.
A second manner to quantify the performance of the adaptive HR approach is by comparing its actual σK and , which were introduced as quality measures in section 3d, with the HR method for each event. The mean values of the σK field () and field () resulting from these techniques are given in Table 3. For the calculation of , only gates with |KDP| ≥ 1° km−1 are considered. Note that from the adaptive HR approach is, on average, one-tenth of the HR method. In addition, results from the adaptive HR approach are smaller than a reasonable error of 20%, while those from the HR method are much larger than 20%. Another quality measure, also given in section 3d, is , which measures the percentage error of the actual μα estimated by the adaptive HR approach. The mean of the field () for each event is found reasonably small as indicated in Table 3. Finally, in terms of computational time required by both techniques, the adaptive HR approach needs, on average, one-third of the time required by the HR method, which is on the order of a few minutes for 1 min of data, while for the conventional technique it is on the order of seconds.
b. During 2 h
Next, the three KDP techniques are compared and evaluated using the same quality measures ρZ,K, , , and and storm events E1–E4 but during 2-h periods as illustrated in Figs. 9–12, respectively. In general, it can be observed that the adaptive HR approach outperforms the other two methods, although the performance of each technique varies according to the storm scenario. For example, the conventional technique can lead to reasonable results when a storm consists of a large area of heavy rain because ΨDP profiles carry sufficient data samples with high SNR levels, reducing the impact from ΨDP outliers. In the HR method, acceptable results can be obtained in a scenario given by a large cell or a cluster of multiple cells with moderate to heavy rain because it allows for consideration of multiple ΔΨDP samples over extended paths, which reduces the impact of small and sometimes negative ΔΨDP values. In addition, this scenario may reduce the sensitivity of the downscaling weight w(i) to noisy measurements of Z and ZDR. In contrast, the adaptive HR approach yields reliable results even when storm cells cover relatively small areas with light rain because ΔΨDP samples are adaptively selected over paths with lengths determined from a predefined interval so that a theoretical value of σK is minimized. Moreover, uncertainties associated with the SC relation are reduced as a result of the improved formulation of w(i) in this method. In case a storm cell becomes significantly small, such that the extents of ΨDP profiles are on the order of Lmin, estimates of KDP by the adaptive HR approach are not possible; this feature could be beneficial because accurate estimation of KDP from limited data samples is rarely achieved.
The resulting time series of ρZ,K for event E3 indicate that during the first hour, the three KDP approaches performed in a similar manner because KDP estimates were obtained from a large cell with heavy rain in which the conventional and HR methods perform at their best. During the second hour, small cells with moderate rain were observed, leading to decreased performance of the conventional and HR methods. In event E2, during 2220–2300 UTC, the HR and adaptive HR approaches provided similar results and performed better than the conventional technique because this period was associated with a cluster of cells with moderate rain. After this period the adaptive HR approach maintained a satisfactory performance, while the performance of the conventional and HR methods decreased because values of KDP were estimated from light rain. For events E1 and E4, ρZ,K time series obtained from the conventional and HR techniques are similar but smaller than those from the adaptive HR approach because in these events single cells with irregular shapes and light rain were observed. Although ρZ,K time series from the three KDP techniques can show similar results for a given storm scenario, Z–KDP relations might include multiple scatters as illustrated in Fig. 8, especially for weak Z values. As a result, the adaptive HR approach provides the best Z–KDP consistency for storm scenarios with different cell sizes and rain amounts.
Time series of from the HR and adaptive HR techniques for the four events exhibited values between 1 and 3 and near 0 (° km−1). However, the results from the HR technique for event E2 showed values smaller than 1° km−1 because of a widespread area of rain with low variability on Z, ZDR, and ΨDP fields. Nonetheless, its performance measured by ρZ,K remains below the adaptive HR approach. This shows a persistent accuracy in estimating KDP by the adaptive HR approach. Furthermore, time series of and resulting from the adaptive HR approach depicted, in all events, consistent percentage errors smaller than 20%. However, for events E1 (around 2200 UTC) and E3 (during some periods after 2040 UTC), the percentage errors increased because of inaccurate measurements of Z and ZDR resulting from storm cells with heavy rain located adjacent to or on top of the radar, leading to power saturation in the receiver. This is an example of how or and or can be used to identify areas where KDP estimates could be compromised. Another example of large can be seen in E2 around 2210 UTC, when values are as high as 40% because of small areas of light rain with a reduced number of ΔΨDP samples, affecting the accuracy of KDP estimates. The discontinuity seen between 2300 and 2320 UTC is due to the constraint |KDP| ≥ 1° km−1. A similar discontinuity is observed in event E1 around 2130 UTC. During the same event, a decreasing and discontinued behavior of ρZ,K is observed in the period 2120–2140 UTC. Such behavior is associated with a progressive reduction of storm cells in intensity and size, which led to light rain echoes with areas smaller than 5 km2 where the extents of ΨDP profiles are not sufficient for the estimation of KDP. This means that KDP cannot be estimated over a ΨDP segment whose length is on the order of or smaller than Lmin. For the case of the HR method, the time series of indicated a limited performance, as values are mostly larger than 50%.
Polarimetric studies have continuously focused on the estimation of KDP because of its capability to overcome power attenuation and radar miscalibration. However, accurate estimation of KDP at scales of the range resolution is challenging because KDP requires significant spatial smoothing because of noisy ΨDP profiles, for example, in light rain. This problem is intensified at short wavelengths when ΨDP profiles include δhv components, for example, in moderate and heavy rain. In this work an adaptive HR approach has been presented to address these problems. The standard deviation of the proposed KDP estimator has been derived and formulated in order to provide a pathlength interval that could lead to KDP estimates with reduced error. This formulation takes into account the spatial variability of the storm and the radar range resolution.
A storm event observed by a polarimetric X-band weather radar during 1 min was used to analyze and test the performance of the KDP estimator. Results showed that the estimated KDP field kept the structure of the attenuation-corrected Z field without significant spatial distortion and that its estimation was associated with reduced errors indicated by the actual standard deviation (i.e., the σK field). The consistency between Z and KDP showed that negative values of KDP, associated with weak Z, can be reduced and that high KDP values, associated with strong Z, can be obtained. To assess the performance of the adaptive HR approach to obtain accurate KDP at range resolution scales, four storm events observed by the same radar but during 2-h periods were considered and the KDP results were compared with the conventional and HR techniques. In general, the proposed approach was able to provide a correlation coefficient between Z and KDP higher than the other two methods. In terms of standard deviations, the adaptive HR approach showed significant improvements compared to the HR technique. The actual μα field, introduced by the adaptive HR approach, was associated with reduced uncertainty as indicated by the results. However, it was observed that the adaptive HR approach is able only to estimate KDP over ΨDP segments larger than Lmin and where the number of ΔΨDP samples is larger than 0.
Although the adaptive HR approach considers measurements of Z and ZDR and constant coefficients related to the SC relation and attenuation correction in rain, the results of this method did not highlight issues related to radar miscalibration, radial patterns in ZDR as a result of PBB, power attenuation, and variability on DSD and drop shape. Consequently, the adaptive HR approach, which uses the correlation between δhv and ZDR and the SC relation, is able to filter δhv and maintain the spatial variability of ΨDP. Therefore, accurate KDP profiles at high spatial resolution in light and heavy rain are achieved. However, if measurements of Z and ZDR are associated with low SNR levels or are affected by offsets that fluctuate along a given radial, then it is expected that the accuracy of the downscaling weights and thereby KDP estimates will be reduced. In general, quality control variables or , which are associated with percentage errors larger than, for example, 20%, might lead to inaccurate KDP estimates.
To achieve the ambitions of implementing the proposed KDP algorithm for real-time operation, further studies are required. This effort includes estimating and updating the coefficients, which are used in the attenuation correction method and in the SC principle, to operational C- and S-band frequencies and, if possible, taking into account temperature variability. Note that at S-band frequency, the ΔΨDP filter condition can be excluded from the KDP algorithm because δhv is usually negligible. Moreover, an automatic algorithm might be needed to classify areas of rain, nonrain, and nonhydrometeors. For solid or mixed rain areas, either the SC ratio should be set to 1, so that KDP is given at expenses of coarse spatial resolution, or another KDP algorithm should be employed. Furthermore, [Lmin; Lmax] should be selected according to values of Δr and to predefined values of μα and σP; the latter two can be set as a function of the rain scenario (e.g., 4 and 3° for convective and 1 and 2° for stratiform rain). Last, it is recommended to evaluate the proposed approach at longer ranges, where attenuation and nonuniform beamfilling affect the Z and ZDR measurements, and in scenarios where PBB is associated with complex terrain features.
A reliable KDP is one of the most powerful observables from polarimetric weather radars. The adaptive HR approach may prove to be key in addressing the dilemma between the spatial resolution and the accuracy of KDP in rain. Moreover, the formulation of the theoretical and the capability of calculating the uncertainty of KDP estimates gate by gate allow for reducing and controlling the errors in the estimation of KDP. Even though this approach still needs to be tested in operational environments, urban hydrology and weather forecast communities may benefit from the proposed approach in terms of spatial resolution, accuracy, and quality control of KDP estimates, which can lead to significant improvements in KDP -based products.
We are grateful to 4TU Centre for Research Data for its support in keeping IDRA data an open access dataset (Russchenberg et al. 2010). The authors extend their deep gratitude to the anonymous reviewers for their constructive reviews, which helped to improve this paper. This work was supported by RainGain through INTERREG IVB NWE Programme.
Standard Deviation of the Estimator
Note that this expression includes the pathlength L and omits the index i for simplicity. Both and α(j) are considered to be independent random variables (rvs) and hence their product is also an rv and denoted by κ(j) with j = 1, 2, … M. Assume that κ(1), κ(2), …, κ(M) have the same variance and that the pair [κ(j), κ(m)], with j ≠ m, has a constant covariance γκ. Using the variance property of the sum of correlated rvs and the relation γκ = , where ρκ is the correlation coefficient, the variance of the KDP estimator is expressed as
If the term is rewritten using and α(j), we obtain
where μα and μΔ represent the mean of α(j) and , respectively. Similarly, and indicate their variance. Equation (A3) can be reduced if ρκ and are assumed to be significantly small. In consequence, Eq. (A3) is simplified to
To include residuals from the ΔΨDP filter condition (i.e., neglecting Δδhv), the difference in differential phase is represented by the sum of two uncorrelated rvs: , where has a mean and variance equal to 0° and , respectively. Thus, Eq. (A4) is rewritten as
For a given path , is expressed as ΨDP(b) − ΨDP(a), which is the difference between two rvs. If both rvs have the same variance , then the variance of is given by = . As a result, a theoretical σK for the KDP estimator is defined as
Filter Design to Estimate by the Conventional Technique
A finite impulse response (FIR) filter type is selected for designing the light filter (first step). The filter order is set to 36, while the required rc is set to 1 km. In addition, a Hann window function is included to obtain a magnitude response H (dB) with small sidelobes. The response H is shown in Fig. B1(a) as a function of the normalized range scale fn = Δr/rs, where rs ≥ Δr and rs in km. Note that H reaches approximately −40 dB at fn = 0.1. This means that the light filter is designed such that spatial variations at range scales smaller than 0.3 km are suppressed. Next, this filter is iterated several times (second step) by setting τ equal to 1.5σp. The storm event E3 at 1950 UTC is used for demonstration purposes. The Z–KDP histograms from the first and second steps are shown in Fig. B1(b). Observe that the iterative step eliminates several outliers without excessively smoothing KDP.
This article has a companion article which can be found at http://journals.ametsoc.org/doi/abs/10.1175/JTECH-D-17-0219.1