## Abstract

Radar-rainfall estimation is a complex process that involves several error sources, some of which are related to the environmental context. The presence of orographic obstacles heavily affects the quality of the retrieved radar products. In relatively flat terrain conditions, dual-polarization capability has been proven either to increase the data quality or to improve the rainfall estimate. The potential benefit of using polarimetric techniques for precipitation retrieval is evaluated here using data coming from two radar systems operating in Italy under complex-orography conditions. The analysis outlines encouraging results that might open new scenarios for operational applications. Indeed, the applied rainfall algorithm employing specific differential phase mostly outperformed the examined reflectivity-based retrieval techniques except for the analyzed winter storm. In the latter case, the likely contamination by frozen or melting snow tended to degrade the performance of the examined *K*_{dp}-based rainfall algorithms.

## 1. Introduction

Rain-rate fields represent valuable information not only for hydrogeological applications but also for microwave communication networks and for assimilation purposes within numerical weather forecast models (Rinehart and Garvey 1978; Joss and Lee 1995; Germann and Joss 2002). In the presence of complex orography, characterized by hilly and mountainous scenarios, rain-rate estimation is a fairly involved process, especially if needed at a ground resolution of less than a few kilometers (Marzano et al. 2004). Rain gauge networks have many limitations related to their sparse and spotlike data distribution. Nevertheless, they represent indispensable means for remote sensor adjustment (often referred to as external calibration) and validation (Vulpiani et al. 2009). Microwave weather radars are considered to be a fairly established technique for retrieving rain-rate fields over large areas from measured reflectivity volumes.

After decades of scientific research following the first pioneering work (Seliga and Bringi 1976), weather radar polarimetry is increasing in popularity within operational meteorological services. Several recent studies and field campaigns have demonstrated the advantages of an operational radar upgrade to polarimetric capabilities. Benefits of weather radar polarimetry include improved data quality, self-consistency techniques for assessment of system calibration (Gorgucci et al. 1999; Gourley et al. 2009), and automatic echo classification (Straka et al. 2000; Marzano et al. 2007, 2008) for weather and forecasting applications.

At frequencies higher than S band, attenuation is the main electromagnetic phenomenon conditioning quantitative precipitation estimation, especially in convective precipitation regimes. Resonance scattering effects (negligible at S band) might also relevantly affect C-band radar measurements (i.e., differential phase or differential reflectivity) and the derived rain rates, even increasing attenuation effects. As a consequence, most experimental studies on the use of dual polarization for operational quantitative precipitation estimation have been performed at S band in relatively flat topographical environments (Ryzhkov et al. 2005b; Giangrande and Ryzhkov 2008; Vulpiani et al. 2009).

A renewed interest in upgrading operational C-band systems to dual polarization showed up since stable polarimetric methods to correct for the effects of attenuation in rain have been introduced through the use of differential phase measurements (Testud et al. 2000; Carey et al. 2000; Bringi et al. 2001; Vulpiani et al. 2008). Most of the European operational weather radars operate at C band—several of them in mountainous areas. In such conditions, together with the enhancement of ground-clutter effects, the major limitation is represented by partial or total beam blocking caused by natural obstructions that very often impose antenna elevations greater than 1.5°. These range-related limitations tend to reduce the potential role of operational weather radars in monitoring precipitation amount at ground level within mountainous areas since, if the nature or intensity of rainfall varies with height (e.g., melting effects during stratiform rain), radar backscattered signals at higher altitudes may not be representative of surface rain rate (Joss and Lee 1995; Germann et al. 2006; Marzano et al. 2004). Moreover, the quality and potential applications of polarimetric radar observations are also sensitive to topography (Friedrich et al. 2009).

The first investigations concerning orographical effects on polarimetric rainfall measurements have been accomplished by Zrnić and Ryzhkov (1996) and Vivekanandan et al. (1999) using S-band radar observations of single flash-flood events. Friedrich et al. (2007) have analyzed the effects of beam shielding on rainfall retrieval at C band in relatively flat topographical conditions for which the main beam blockage source was represented by urban obstacles. Silvestro et al. (2009) have more recently proposed and analyzed the performance of a flowchart–decision-tree polarimetric algorithm applied to radar measurements collected with an operational C-band system located at 1400 m above mean sea level (MSL) in the western Apennine Mountains (Liguria region, Italy). The latter study, restricted only to unshielded sectors, to 100 km of range, and to daily rainfall evaluation, has shown that polarimetric variables, whenever used by the algorithm (only about 10% of the total cases), improved rainfall estimation.

The work presented here is aimed at investigating, without any constraints (i.e., range distance or shielded sectors), the potential benefit derived by the use of dual polarization for operational precipitation estimation in complex-orography scenarios. For this purpose, we employed data acquired by the C-band dual-polarized radars located at Il Monte mountain near Tufillo (Chieti, Abruzzo region, central Italy) and Zoufplan mountain near Paluzza (Udine, Friuli Venezia Giulia region, northern Italy) and by dense networks of more than 150 rain gauges within each radar’s nominal coverage of 175 km. Ten days of observations, for a total of 240 h, were analyzed. In particular, the work is aimed at 1) clutter removal through the synergy of clutter map, radial velocity, and polarimetric texture analysis; 2) correction of two-way path attenuation; 3) reconstruction of the vertical profile of reflectivity; and 4) polarimetric estimation of near-surface rainfall. Several techniques to address these problems are intercompared and evaluated. The paper is organized as follows. Section 2 provides information on the employed dataset. The data processing chain is described in section 3. The results are discussed in section 4. Details of the proposed technique for differential phase processing are provided in the appendix.

## 2. Dataset description

### a. Radar systems

The radar systems used in this work belong to the Italian network that provides coverage of most of the Italian territory with updates every 15 min and ground spatial resolution lower than 1 km. The primary justification for a weather radar network in Italy, as well as in other countries, is the detection and warning of severe weather and related hydrogeological risks. Hydrological risk is further enhanced by orography, which is characterized in Italy by small catchments along most coastlines and by the Alpine and Apennine chains. Because of the incomplete national coverage, the Department of Civil Protection (DPC) has been appointed to be responsible for complementing and integrating the existing systems. These systems are made of 10 C-band installations belonging to regional authorities, five of which are polarimetric, and one transportable X-band polarimetric radar. Two systems are owned by the Italian company for air navigation services (ENAV), and three are managed by the Meteorological Department of the Italian Air Force (AMI). The Polarimetric Doppler Radar Systems (PDRS), located at Il Monte (hereinafter PDRS1) and Zoufplan (hereinafter PDRS2), respectively, are two of the six DPC polarimetric radars. PDRS1 and PDRS2 are respectively located near the border between the Molise and Abruzzo regions and in the Friuli Venezia Giulia region, close to the border with Austria, as shown in Fig. 1.

The Italian radar network siting resulted from a compromise between the radar network fulfilling needs and the logistics and environmental requirements. The PDRS1 location at about 700 m MSL is surrounded on the northwestern side by the highest peak (about 3000 m) of the Apennine mountain range (see Fig. 1). The main beam blocking toward the inland country is caused by the Maiella mountain (the highest peak is at about 2800 m) at a distance of about 35 km. In such circumstances, the major error sources in radar-rainfall estimation are obviously related to ground-clutter contamination, partial and/or total beam shielding, and vertical variability of reflectivity. The PDRS2 site (at about 2000-m altitude) has been chosen to observe the Friuli Venezia Giulia and Veneto valleys, which are characterized by the presence of several catchments (e.g., Tagliamento, Livenza, and Isonzo) that are almost fully visible even at low antenna elevation scans.

The considered PDRSs operate at 5.6 GHz with a maximum unambiguous range of approximately 175 km and a range resolution of 0.15 km. The antenna rotation speed has been set at 12° s^{−1}, allowing an integration of 68 samples within the resolution volume. The two systems are connected by satellite links to the two National Radar Primary Centres (RPC), located in Rome at DPC and Savona at the Centro Internazionale in Monitoraggio Ambientale (CIMA) Research Foundation to mainly ensure remote control, through the Radar Remote Control server, and products generation. The RPC at Savona works as a “backup center” to guarantee that the system is continuously operating. The Radar Archive Centre subsystem is devoted to archiving and managing radar data and products by means of a relational database. The generated products are then disseminated to all institutions composing the national network. The considered radar systems are subject to periodical in situ maintenance activities (at least four times per year) during which the calibrations of the receiver and transmitter are checked. In addition, the receiver calibration is also remotely performed through specific maintenance software.

### b. Test cases

Seven precipitation events for a total of 10 days of observations have been analyzed in this study. The events observed by PDRS1 and PDRS2 are respectively listed in Table 1 and Table 2 together with the maximum observed cumulated rainfall over different time periods, the number of available gauges, and the freezing-level height (FLH) as derived from the closest available radiosoundings.

Most of the considered events occurred between spring and autumn, and the 2-day event observed by PDRS1 at the beginning of March of 2011 is a typical winter storm of the Mediterranean Sea region. Indeed, the Euro-Atlantic synoptic configuration favored cold polar fluxes into the Mediterranean area where a deep cutoff low was acting. Specifically, the center of Italy was affected by a cold front with an FLH that varied between 1 and 1.5 km MSL.

The 3-day weather event that occurred at the beginning of May of 2010 (observed by PDRS2) was the most interesting in terms of persistence and severity. It was characterized by a deep low pressure area (i.e., geopotential height of about 5442–5460 m at 500 hPa) over the central-western Mediterranean region, as can be deduced from Fig. 2, which shows the superposition of the geopotential height at 500 hPa and the mean sea level pressure [from European Centre for Medium-Range Weather Forecasts (ECMWF) numerical model analysis]. The Euro-Atlantic zone was characterized by a very low zonal index with an Ω-blocking structure. At the beginning of the event the Italian territory was affected by a warm conveyor belt transporting wet warm fluxes from the African continent. The crucial phase of the storm happened on 4 May with the beginning of the cyclogenesis on the Gulf of Lion (between France and Spain, close to Marseille), as depicted in Fig. 2b.

The blocking structure caused the persistence of precipitation with convective phenomena in central-northern Italy, especially over the Alps and northern Apennine areas. The cumulated rainfall field derived from interpolated gauge observations is shown through contour lines in Fig. 3. The peaks of cumulated precipitation have been observed over mountain areas, within about 60–70 km from the radar site. The maximum 72-h cumulated rainfall was about 270 mm, with a maximum daily accumulation of about 153 mm on 4 May 2010 (see Table 2 for details).

## 3. Data processing

This section describes the radar-data processing chain applied to the observed precipitation events (listed and described in section 2b). As depicted in Fig. 4, the processing begins with the identification of and compensation for nonmeteorological echoes and partial beam blockage effects. Second, *K*_{dp} is derived from the filtered Φ_{dp} and used in the rain path attenuation correction module, where external data (i.e., temperature soundings) are also ingested. The mean vertical profile of reflectivity (VPR) is then retrieved and applied to the considered reflectivity-based radar products [i.e., the lowest beam map (LBM) and vertical maximum intensity (VMI)] before computing the rainfall rate. The latter is also estimated through the retrieved specific differential phase.

### a. Artifact removal and partial beam-blocking correction

As mentioned in section 2a, both PDRS are surrounded by high mountains; as a consequence, the major error sources in radar-rainfall estimation are obviously related to ground-clutter contamination, partial and/or total beam shielding (expecially for PDRS1), and altitude of the measurement above ground (i.e., vertical variability of reflectivity).

The radar echoes generated by nonmeteorological targets are discriminated from weather returns by means of a quality map *Q* that is subjectively generated by combining the following quality indicators: static-clutter map (CMAP), radial velocity *V*, texture of *Z*_{dr} (TxZdr), *ρ _{hυ}* (TxRho), and Φ

_{dp}(TxPhi). CMAP is a volumetric map obtained by averaging a wide set of reflectivity data observed in clear-air conditions. For each quality indicator

*X*(i.e.,

_{j}*X*

_{1}= CMAP,

*X*

_{2}=

*V*,

*X*

_{3}= TxZdr,

*X*

_{4}= TxRho, and

*X*

_{5}= TxPhi), the degree of membership in the nonmeteorological target class is defined through a trapezoidal transformation function

*d*=

_{j}*d*(

*X*):

_{j}where *X _{i}*

_{,j}is the

*i*th vertex of the trapezoid relative to the

*j*th quality indicator. Table 3 shows the parameterization used for defining

*d*. The relative quality index

_{j}*q*associated with

_{j}*X*is then defined as the complement of

_{j}*d*(i.e.,

_{j}*q*= 1 −

_{j}*d*). Last, the overall quality

_{j}*Q*is obtained through a weighted sum of the relative quality indices:

Radar returns with associated low quality (i.e., *Q* < 0.5) are finally rejected.

To take into account the beam-shielding effects properly, an electromagnetic propagation model is used to identify the obstructed radial directions. The partial beam blockage (PBB) map, representing the occultation degree at a specific antenna elevation, has been retrieved by resorting to the simplified obstruction function proposed by Bech et al. (2003), assuming the wave propagation in the standard atmosphere (Doviak and Zrnić 1993). The estimated PBB is then compensated up to 70%, as in Tabary (2007).

Figure 5 shows the visibility maps (i.e., the complement of the PBB maps) in terms of percentage of nonoccluded beam for PDRS1 (top panel) and PDRS2 (bottom panel). As can be seen in the top panel of Fig. 5, the visibility of PDRS1 is less than 200° wide at 0.4° of antenna elevation. This situation is particularly critical considering that only a reduced coast area is visible at low-elevation scans, whereas the internal region becomes fully visible only at about 2.5°. Considering that PDRS2 is sited at about 2000 m of altitude, the visibility of PDRS2 is basically guaranteed southward, where the main settlements are located. The partial or total northward shielding (the Alps area) does not represent a concern.

### b. Differential phase processing

As described in the appendix, a polarimetric radar system provides measurements of the total differential phase Ψ_{dp} that is the sum of the differential propagation Φ_{dp} and the backscatter phase *δ _{hυ}*. We are only interested in the propagation component for attenuation correction and rainfall estimation purposes,

*K*

_{dp}being related to the range derivative of Φ

_{dp}. At C-band frequencies,

*δ*might not be negligible when resonance scattering occurs. Furthermore, Ψ

_{hυ}_{dp}is also conditioned by system noise, offset, and potential aliasing problems. Several approaches have been proposed to handle this issue, for example, median filter, moving average, polynomial fitting (Bringi and Chandrasekar 2001), finite impulse response filter (Hubbert and Bringi 1995), and more recently range derivative in the complex domain (Wang and Chandrasekar 2009). In the work presented here, a multistep moving-window range derivative approach is applied. As described in Fig. 6, the applied method can be summarized in four main steps:

*K*_{dp}check—special care is taken in treating the*K*_{dp}values that are not manifestly physical,*K*_{dp}typically ranging between −2° km^{−1}(as for vertically oriented ice crystals) and 20° km^{−1}(as for heavy rain) at C band (Bringi and Chandrasekar 2001);*K*_{dp}retrieval (final guess)—the final estimation of the specific differential phase is then obtained from .

In step 2, we attempt to discriminate anomalous *K*_{dp} values coming from aliasing or other phenomena (i.e., noise, backscatter differential phase, nonuniform beam filling, or residual artifacts). For the latter cases the estimated *K*_{dp} are set to zero, but in the case of Φ_{dp} wrapping it is necessary to refine the retrieval process. Otherwise, according to step 3, the reconstructed profiles of differential phase shift would contain *L*-sized range segments characterized by unrealistic stationary Φ_{dp} values. Indeed, when aliasing occurs, Φ_{dp} is exposed to a folding of the same order of magnitude as the maximum unambiguous phase shift Φ_{dp,max} (i.e., 360° for simultaneous transmission). As a consequence, the estimated *K*_{dp} values would be systematically low (i.e., *K*_{dp} ~ −0.5 × 360/*L* ~ −25° km^{−1}, with *L* = 7 km) in any *L*-sized range segment centered at range *r _{a}*, where aliasing shows up.

Once any of such range segments is identified, Φ_{dp} is unwrapped by adding Φ_{dp,max} and the whole processing procedure is repeated from *r _{a}* −

*L*/2 to the end of the range profile. It is interesting to note that following steps 1–4 the retrieved differential phase is not affected by the system offset, it being removed when computing .

It is also worth mentioning that, as shown in the appendix, the standard deviation of becomes

where *N* is the number of range gates contained in the *L*-sized moving window (i.e., *N* = *L*/Δ*r*, Δ*r* being the range resolution). This means that, for *L* = 7 km and Δ*r* = 150 m and assuming *σ*(Ψ_{dp}) = 3°, is about 0.05° km^{−1}.

As outlined in the appendix, the general expression for the standard deviation of *K*_{dp} estimated by recursively iterating steps 3 and 4 is

where *I* is the number of iterations (with *I* ≥ 1). Consequently, for a lower window size and/or a poorer range resolution Δ*r* (i.e., lower *N*), it might be possible to get about the same standard deviation by just iterating the retrieval procedure a few times [e.g., for *L* = 4 km and Δ*r* = 300 m with *I* = 2, assuming *σ*(Ψ_{dp}) = 3°].

As an example, the top panel in Fig. 7 shows the range profile of the raw (dashed line) and filtered differential (dotted line) phase. To evaluate better the effect of the filtering approach, the current system offset has also been added to the filtered Φ_{dp} (see the continuous curve). The bottom panel in Fig. 7 shows the corresponding retrieved specific differential phase. A direct assessment of the quality of the applied *K*_{dp} retrieval approach is beyond the scope of this work. The quality of the specific differential phase reflects on the corresponding rainfall rates that are evaluated in depth later in the paper, however. Nevertheless, the effectiveness of any *K*_{dp} estimation technique can only be evaluated qualitatively by resorting to the consistency between estimated and theoretical values of *K*_{dp} (as obtained from scattering simulations) in a *K*_{dp}-versus-*Z _{hh}* plane, as was done in Wang and Chandrasekar (2009).

### c. Attenuation correction

Rain path attenuation is known to be one of the main impairments in radar-rainfall estimation at frequencies higher than S band. Stable procedures for accounting for that issue are available through the use of differential phase shift. Indeed, specific differential phase and specific attenuation are almost linearly related (Bringi et al. 1990):

where *α _{hh}*

_{,dp}is the copolar specific attenuation and specific differential attenuation (dB km

^{−1}), respectively. Using (4), the two-way path-integrated attenuation

*A*

_{hh}_{,dp}can be written as

and the corrected reflectivity and differential reflectivity are

where *r*_{0} identifies the beginning of the rain path. Despite the simplicity of (5)–(6), it is necessary to deal with the space–time variability of *γ _{hh}*

_{,dp}, they being related to raindrop temperature, size, and shape (Jameson 1992). Carey et al. (2000) suggested retrieving the optimal values of

*γ*

_{hh}_{,dp}by analyzing the trend of the observed

*Z*

_{hh}_{,dr}with respect to Φ

_{dp}after minimizing their intrinsic variability. In Vulpiani et al. (2008) a preliminary rain-type (or, in other terms, average raindrop size distribution) identification enables one to estimate

*γ*

_{hh}_{,dp}on a gate-by-gate basis. In place of the Bayesian approach used in Vulpiani et al. (2008), a fuzzy-logic hydrometeor classification algorithm, adapted from Marzano et al. (2007) by also including the correlation coefficient as input, has been embedded within the adaptive optimization procedure. As may be expected, the accuracy of the classification algorithm degrades in the presence of any bias in the input variables. Marzano et al. (2007) estimated an accuracy of about 70% in the presence of a bias of 2.5 dB on

*Z*and 0.75 dB on

_{hh}*Z*

_{dr}. Therefore, it is reasonable to expect that, with a standard 1-dB uncertainty on the reflectivity (the system is reasonably well maintained) and also making use of the correlation coefficient, a bias on differential reflectivity, estimated to not exceed a few tenths of decibels for both PDRS1 and PDRS2, should not significantly impair the algorithm performance.

As depicted in Fig. 8, the overall adaptive Φ_{dp} correction procedure can be summarized through the following few steps:

The radar observables are filtered from nonmeteorological targets and the phase measurements are processed as described in section 3b.

A preliminary attenuation correction is performed using (5)–(6) assuming fixed values for

*γ*_{hh}_{,dp}(i.e., 0.08 and 0.02 dB km^{−1}). At this stage the temperature profile*T*, retrieved from the closest available radiosounding, is used to roughly discriminate rain from frozen particles.The corrected

*Z*_{hh}_{,dr}are then used with*K*_{dp},*ρ*, and_{hυ}*T*for hydrometeor classification.Values of

*γ*_{hh}_{,dp}are associated with each rain type (i.e., light, moderate, heavy, or large drops) as derived from scattering simulations (Vulpiani et al. 2008).

### d. Rainfall estimation

Rainfall rate is here computed from reflectivity and specific differential phase. It is known that *Z*_{dr} is another potential rain-rate predictor when coupled with *Z _{hh}* or

*K*

_{dp}(Bringi and Chandrasekar 2001). The observed

*Z*

_{dr}has been found to be azimuthally biased by nearby obstacles (i.e., lightning rods, or fencing), however. Furthermore, any

*Z*

_{dr}-based rainfall estimation algorithm is generally more sensitive to uncompensated rain path and wet radome attenuation effects.

Once reflectivity is corrected for attenuation, a mean VPR (Joss and Lee 1995) is retrieved from every volume scan (provided that it contains meteorological echoes at all defined height levels) and is then applied either to the LBM or to the VMI to get ground-projected reflectivity products. We have also considered nonprojected reflectivity fields as reference radar products. A standard *Z*–*R* relationship (Marshall and Palmer 1948) is applied to reflectivity products. For the *K*_{dp}-based rain-rate algorithm, we have considered a general expression of the form *R* = *a*|*K*_{dp}|* ^{b}* sign(

*K*

_{dp}), as suggested by Ryzhkov et al. (2005a). Although this formula provides unrealistic negative rain rates for negative

*K*

_{dp}values, it is adopted to compensate for the noise effects on the retrieved

*K*

_{dp}(i.e., slightly positive and negative rain rates tend to cancel each other when computing the cumulated rainfall). Considering the power-law parameters

*a*and

*b*derived by either Bringi and Chandrasekar (2001) or Scarchilli et al. (1993), the tested algorithms, denoted as

*R*

_{BC}and

*R*

_{SC}, respectively, are

where *f* is the radar frequency expressed in gigahertz. Relationships (7) and (8) have been applied to the lowest beam map of *K*_{dp} (LBMK). Considering the characteristics of the proposed *K*_{dp} estimation technique [i.e., *σ*(*K*_{dp}) ≈ 0.05° km^{−1}], it makes sense to compensate rain rates corresponding to −0.05° km^{−1} ≤ *K*_{dp} ≤ 0.05° km^{−1}. Wherever *K*_{dp} drops below −0.05° km^{−1}, the average of the neighboring values (within an area of 9 km^{2}) exceeding such a threshold is considered in its place. Otherwise, if none of the neighboring values satisfies such conditions, the rainfall rate is set to zero.

The following notation is used to identify the considered algorithms:

*R*[LBM(*Z*)] for rainfall rate computed from LBM,*R*(LBM_{VPR}) for rainfall rate computed from ground-projected LBM (through the mean VPR),*R*[VMI(*Z*)] for rainfall rate computed from VMI,*R*(VMI_{VPR}) for rainfall rate computed from ground-projected VMI (through the mean VPR),*R*_{BC}(LBMK) for rainfall rate computed from LBMK using (7), and*R*_{SC}(LBMK) for rainfall rate computed from LBMK using (8).

## 4. Results

As usual, the performance analysis of the considered rainfall algorithms is accomplished by comparing radar-rainfall fields with in situ measurements provided by the available gauge network. Considering the intrinsic difference between the two measurement sources, the radar–gauge matching criterion might not be an irrelevant feature of the validation step, especially in complex-topography scenarios in which the radar volume altitude is typically high. In such circumstances, wind drift might represent an important effect to address. Similar to the approach followed by Silvestro et al. (2009), the radar estimation minimizing the square difference with respect to the gauge observation has been selected within a 5 × 5 km^{2} area around the gauge position. Although this method might lead to consideration of different radar volumes for each evaluated rainfall algorithm, the best estimation selection is always guaranteed for all of them.

This section first describes a qualitative performance analysis to get a general overview; then, a quantitative statistical evaluation of the hourly cumulated precipitation is discussed in terms of the following error indicators:

mean error ,

error standard deviation ,

root-mean-square error, defined as ,

mean bias, defined as the mean ratio between gauge observation

*R*and radar estimate_{G}*R*—that is, bias = 〈_{R}*R*/_{G}*R*〉, and_{R}correlation coefficient (CC), defined as ,

where *R _{R}* is the radar estimate and

*R*denotes the gauge observation. The angle-bracket notation stands for the average operator with respect to space and/or time, and represent the standard deviations of the gauge and radar observations, respectively.

_{G}### a. Overall analysis

Figures 9 and 10 show the spatial average of the cumulated rainfall as a function of time for all of the considered algorithms relative to the storm events observed by PDRS1 and PDRS2, respectively. These figures indirectly provide information on the mean error as a function of the accumulation time. Except for events VI–VII, depicted in Fig. 9f, it can be noticed that *R*_{BC} generally outperforms *R*_{SC} and all of the methods that employ reflectivity when compared with reference gauges (blue curves in Figs. 9 and 10). As long as the cumulated precipitation does not exceed 5–10 mm, *R*(LBM_{VPR}) generally provides a relatively good estimation; for example, cases II and III, respectively shown in Figs. 9b and 9c, are emblematic.

Otherwise, both of the considered *K*_{dp}-based rainfall algorithms depart from the observed rain depths by more than *R*(VMI_{VPR}) and *R*(LBM_{VPR}) do for the analyzed winter event (cases VI and VII). Indeed, *R*_{BC,SC}(LMBK) are likely conditioned by frozen hydrometeors while the reflectivity-based techniques, especially *R*(LBM_{VPR}), take benefit from the ground projection.

As can be noticed from Fig. 10, *R*_{BC}(LBMK) provides, on average, the best performance for the events observed by PDRS2, whereas *R*_{SC}(LBMK) and especially the reflectivity-based techniques clearly underestimate precipitation.

Figure 11 shows the spatial distribution of the mean bias computed on the hourly cumulated rainfall maps for test cases I, VI–VII, and VIII–X relative to *R*[VMI(*Z*)], *R*(VMI_{VPR}), *R*(LBM_{VPR}), and *R*_{BC}(LBMK). In general, *R*_{BC}(LBMK) provides relatively uniform bias fields with a tendency to underestimate precipitation in the shielded sector mainly at longer ranges because of the likely contamination by frozen particles. For the winter events (VI–VII), it can also be noticed that *R*_{BC}(LBMK) slightly overestimates precipitation in the visible sector, in agreement with the findings shown in Fig. 9. Short-range underestimation is also evident in the same case, however. This result might be attributed to a disturbance in the measured differential phase caused by sidelobe contamination of low-intensity precipitation observations. As can be seen in Fig. 11, the *R*[VMI(*Z*)] algorithm generally underestimates rainfall in the blocked sectors where VMI(*Z*) is constructed from higher elevation scans that might cause precipitation overshooting and/or contamination by ice particles.

The VMI(*Z*) ground projection through the application of the mean VPR enables one to mitigate these effects, however. The improvement determined by the VPR correction is outstanding for all considered events, especially for the analyzed winter storm relative to *R*(LBM_{VPR}). Moreover, the VPR correction causes an average overestimation on the VMI-based rainfall algorithm for cases II and III, as shown in Figs. 9b and 9c, respectively.

### b. Error analysis

The overall results discussed in section 4a are quantitatively confirmed through the considered error statistics summarized in Table 4 for *R*(VMI_{VPR}) and *R*_{BC}(LBMK) relative to PDRS1. Except for cases I, IV, and V, the performance of *R*(VMI_{VPR}) with respect to *R*_{BC}(LBMK) is characterized by a comparable mean error and mean bias, and the correlation is generally lower for *R*(VMI_{VPR}). It is unmistakable that *R*_{BC}(LBMK) outperforms *R*(VMI_{VPR}) in terms of RMSE by virtue of a lower error standard deviation. This behavior, which is particularly evident for events II, IV, and VI–VII, might be attributed to the ground projection by means of the mean profile of reflectivity, which is unable (by construction) to capture the VPR spatial variability. Of interest is that *R*_{BC}(LBMK) produces a mean bias that is very close to the optimal value for test cases I and V, with relatively small deviations from unity for the remaining events. In addition, it is important to mention that *R*_{BC}(LBMK) seems to work efficiently even for a moderate storm such as that observed on 23 October 2009 (case V).

To evaluate the beam-blocking effects on rainfall retrieval, the error analysis has also been limited to the shielded and unshielded azimuthal sectors. As can be noticed by comparing Tables 5 and 6, the beam blockage and the corresponding need to use high elevation scans are the main error source for the considered reflectivity-based rainfall algorithms. Otherwise, *R*_{BC}(LBMK) provides better error scores in the shielded areas as compared with the unshielded areas for some of the analyzed cases (i.e., I and VIII–X). This effect might be related either to the storm characteristics or to the intrinsic error structure of the applied rainfall technique. Indeed, a general tendency to overestimate precipitation, as can be deduced from Fig. 9, might be compensated by the blocking effects themselves.

For the 3-day event observed by PDRS2, the performance of *R*(VMI_{VPR}) is worse than *R*_{BC}(LBMK) especially in terms of *σ _{ɛ}* and mean bias (see Table 7). More specifically,

*R*(VMI

_{VPR}) tends to generally underestimate the precipitation whereas

*R*

_{BC}(LBMK) does not show any systematic trend.

Figure 12 shows the range dependency of the mean bias in the shielded (left panels) and visible (right panels) azimuthal sectors for the rainfall algorithms *R*[VMI(*Z*)] (top panels), *R*(VMI_{VPR}) (middle panels) and *R*_{BC}(LBMK) (bottom panels) relative to PDRS1 and test case I. On one hand, the VMI-based algorithms are characterized by an increasing bias at longer distances in the obstructed sector, although they are generally mitigated by the ground projection. On the other hand, all of the rainfall methods seem to be insensitive to range distance when considering the visible region. This effect might be partially justified given the storm characteristics (i.e., FLH), the altitude of the radar site (i.e., 700 m MSL), and the scan strategy. Indeed, because of the fact that the unshielded sector is visible at a low antenna elevation angle (i.e., 0.4°), the contamination by melting or frozen snow (among the main sources for the vertical variability of reflectivity) was unlikely up to a distance of about 120 km for the analyzed case. The same behavior is confirmed by event V, as depicted in Fig. 13.

The range dependency of the rain fields estimated by *R*[VMI(*Z*)] is more evident for the 3-day event observed by PDRS2, the radar being sited at about 2000-m altitude. Indeed, Fig. 14 outlines that the considered reflectivity-based algorithm takes advantage of the VPR correction even in the unshielded azimuthal sector. As already noted in Fig. 11, *R*_{BC}(LBMK) shows a slight trend toward underestimation in the blocked sector, especially at far ranges.

The overall fair efficiency of *R*_{BC}(LBMK) can also be deduced from Figs. 15 and 16, showing the comparison among gauge and radar estimates of the total cumulated rainfall respectively for PDRS1 and PDRS2. The error statistics are also summarized for each event in terms of RMSE, mean bias, and correlation coefficient. For completeness, the average value of the measured precipitation 〈*R _{G}*〉 is also shown to evaluate quickly the fractional standard error (i.e., FSE = RMSE/〈

*R*〉). As can be noticed, although most of the data are distributed around the bisector, the largest bias (i.e., bias = 1.47) has been found for event IV, characterized by FSE(%) ≈ 40. The highest FSE (i.e., ≈56%) was found for event III, however, which was marked by relatively low mean precipitation. Moreover,

_{G}*R*

_{BC}(LBMK) performed reasonably well (i.e., FSE lower than about 30%) for events I, V, and especially VIII–X.

## 5. Conclusions

The potential benefit derived from the use of polarimetric methods for operational precipitation estimation in complex-orography scenarios has been investigated. A couple of dual-polarized C-band radars (named PDRS1 and PDRS2) belonging to the Italian radar network have been considered together with a dense gauge network. PDRS1 is located in central-eastern Italy near the coastline and close to the Apennine range (at about 700 m MSL), and PDRS2 is sited in the middle of the northeastern Alps (at about 2000 m MSL). The local climate and observation geometry of the two sites are very different, and in this respect they represent typical cases of operational radar in complex orography. Spring and autumn mesoscale precipitation systems are mainly analyzed because those are the rainiest seasons in northern and central Italy.

Most of the error sources affecting operational radar-rainfall estimation have been addressed. A combination of clutter-map, radial-velocity, and polarimetric-texture analysis is applied for the evaluation of data quality to suppress nonmeteorological echoes (i.e., ground clutter, clear-air echoes, and interferences caused by wireless local area networks). Partial beam-blocking effects are accounted for by resorting to an electromagnetic propagation model based on a 240-m digital elevation model. A new fairly efficient algorithm for differential phase measurement processing and specific differential phase estimation is applied. Rain path attenuation effects are also handled through the adaptive use of differential phase measurements. Rain-rate fields are finally retrieved either from reflectivity-based radar products (LBM or VMI, eventually ground projected through the estimate of the mean profile of reflectivity) or from the specific differential phase.

The comparative analysis accomplished for five 1-day events, one 2-day event, and one 3-day event has shown promising and valuable outcomes from the use of *K*_{dp} for operational precipitation estimation in the presence of complex environmental scenarios. Rainfall fields estimated from the lowest beam map of specific differential phase have generally matched the rain gauge observations better, especially in the shielded areas. It is worth mentioning that these results are also confirmed for low-to-moderate rain rates in all considered cases. This may suggest that a *K*_{dp}-based algorithm may be even applied without resorting to a decision tree where *Z _{hh}*-based algorithms are employed for less-intense rainfall (with the open problem of how to decide the geographically and storm-dependent threshold values).

This work has outlined the sensitivity of the considered *K*_{dp}-based rainfall estimators to the presence of dry or melting ice particles. Therefore, future work will analyze the sensitivity of polarimetric radar observables on precipitation regime and seasonal dependency to set up a rain retrieval technique that could adaptively be used for operational purposes.

## Acknowledgments

This work is partially funded by the Italian project IDRA commissioned by the Department of Civil Protection to CETEMPS (University of L’Aquila). The authors are grateful to Prof. B. De Bernardinis, Dr. P. Pagliara, engineers M. Negri, L. Rossi, and F. Santamaria, and Prof. F. Siccardi for the efforts devoted to the realization of the national radar network.

### APPENDIX

#### Accuracy of the Specific Differential Phase

The retrieval process of *K*_{dp} from the measured differential phase Ψ_{dp} might not be an easy task. Indeed, the *K*_{dp} estimation is conditioned by the backscatter differential phase *δ _{hυ}* that adds to the propagative component Φ

_{dp}:

Therefore, the accuracy of the *K*_{dp} retrieval depends either on the approach adopted to derive it from the measured differential phase or on the system noise affecting Ψ_{dp}.

In general, the variance of *K*_{dp} can be estimated by resorting to the error propagation theory (Taylor 1997). In the case in which *K*_{dp} is estimated through a finite-difference scheme on moving window of size *L*,

the resulting variance is (Bringi and Chandrasekar 2001)

where Δ*r* is the range resolution and *N* is the number of range gates contained in the moving window.

As described in section 3b, the proposed algorithm for *K*_{dp} estimation and Φ_{dp} reconstruction consists of four steps. At step 3, the differential phase is obtained as a range integral of the specific differential phase retrieved at step 1 (called ). It clearly means that

Using (A4), the final specific differential phase that is computed as a range derivative of becomes

According to the uncertainty propagation theory, the variance of a linear combination of independent random variables is equal to the linear combination of the corresponding variances. As a consequence, the variance of can be written as

Using (A3), the standard deviation of becomes

Taking into account (A6), it is easy to verify that the standard deviation of the specific differential phase might be further reduced by iterating steps 3 and 4:

where *I* is the number of iterations (with *I* ≥ 1).

It is worth noting that approximating Φ_{dp} through a linear regression would lead to the following expression for the standard deviation of *K*_{dp} (Bringi and Chandrasekar 2001):

which is bigger by a factor of (6)^{1/2} than , as expressed in (A7).