Evaluations on radar QPE using raindrop size distribution in Southern Luzon, Philippines

The study analyzed the raindrop size distribution (DSD) measured by an optical Parsivel disdrometer in Southern Luzon, Philippines and utilized it to generate dual-pol relations for the nearby Tagaytay radar. The relations were generated using two methods (Method 1 - gamma-based and Method 2 - linear fitting), four time-integration steps (1, 2-, 5, and 10-min) and datasets from two periods (wet season and single event). The resulting quantitative precipitation estimates (QPEs) calculated from the generated R(Z) relations were compared to rain gauge stations near the disdrometer and were evaluated for the Tropical Storm Yagi Monsoon event of 10 August (2200 UTC) to 11 August (0400 UTC) 2018 using six statistics: Pearson’s correlation; mean error, percent bias, Nash-Sutcliffe Efficiency, mean absolute error, and root-mean-square error. Results show that the area’s DSD demonstrates relatively larger average raindrop diameters than some of its Asian counterparts, albeit a smaller number in the total number of raindrops when compared with the same areas. In terms of QPE evaluation, results showed a consistent pattern observed wherein the R(Z) relations using finer time steps (1-and 2-min) generally performed better than the longer ones. Moreover, Method 1 dominated Method 2 in terms of error statistics. As expected, Method 2 outperformed Method 1 in terms of r (as Method 2 itself is derived through linear fit). The best derived R(Z) relations were able to outperform other relations in terms of r, NSE, and RMSE. On the other hand, R(K DP ) was able to perform the best in terms of ME, MAE, and pBIAS, reducing the bias of current standard method by up to 74%.


INTRODUCTION
Quantifying precipitation over time constitutes a number of significant limitations due to its highly discrete and variable nature (Krajewski et al. 2003;Kathiravelu et al. 2016).As a major driving force of global water cycle, precipitation data is extremely valuable in hydrology, meteorology, agriculture, and weather forecasting (Trenberth et al. 2003;Kidd and Huffman 2011) and with the continuous progress of technology, the sheer availability of rainfall data is not enough (Jiang et al. 2012).The need for high-quality, comprehensive, and accurate precipitation measurements and estimates is becoming more and more important for hydrological modeling (Liu et al. 2017;Lee et al. 2019a), numerical weather prediction (Rossa et al. 2008;Shrestha et al. 2013), and other hydrometeorological applications (Yilmaz et al. 2005;Pappenberger et al. 2008).
Measuring rainfall is performed using three conventional methods: (1) direct measurement via rain gauge; (2); remote sensing through satellites; and (3) use of groundbased weather radars.Direct measurement of rainfall through rain gauges is widely used and is often considered as true and reference rainfall (Richards and Crozier 1983).Terr. Atmos. Ocean. Sci., Vol. 32, No. 5, Part I, 693-724, October 2021 A major disadvantage, however, is the limited spatial range of a rain gauge which is highly influenced by the surrounding features (Rana et al. 2015;Kidd et al. 2017).The obvious solution to overcome this is the deployment of a dense network of rain gauges that captures rainfall over a specifically homogenous region.However, this is extremely costly especially in the mountainous regions where the density of rain gauges should theoretically be the highest.Furthermore, errors have the propensity to increase in magnitude with respect to the degree of complexity of a terrain and rainfall conditions (Daly et al. 1994).Besides the use of conventional rain gauges, the use of satellite products is also a growing field in generating precipitation data (Skomorowski et al. 2001).Satellite rainfall usually provides precipitation data at a global scale at various spatial resolutions (Xie et al. 2003;McPhee and Margulis 2005;Joshi et al. 2013).However, its use is limited to basins large enough to fit through multiple grids depending on the satellite data's resolution.Moreover, most satellite-derived rainfall datasets are not always available in real-time and can be difficult to access during flash floods where the need for rainfall data is at its peak significance (Grose et al. 2002).The use of Doppler weather radars for measuring rainfall not only provides a spatial range large enough to include several regions at a time but also offers real-time precipitation information in high spatial and temporal resolution (Scofield and Kuligowski 2003;Germann et al. 2006;Yoon et al. 2012;You et al. 2018).Instead of measuring rainfall directly, radars retrieve rainfall data by emitting pulses of electromagnetic energy at microwave frequencies and measuring the resulting reflectivity (Z) scattered by raindrops in the atmosphere within a scanning volume and is usually measured in decibel relative to Z or dBZ.Studies show that weather radars may be preferred over rain gauges for three major reasons (Richards and Crozier 1983;Vieux and Bedient 1998;Löffler-Mang and Blahak 2001;Chandrasekar and Cifelli 2012;Dutta et al. 2012): (1) it can include remote and inaccessible regions where the installation of rain gauges may prove nearly impossible; (2) it can cover areas where the density of rain gauge is insufficient to characterize the region; and (3) the provision of 3-dimensional areal coverage eliminates the need for extrapolating point measurement over an area of interest.
In the Philippines, a network of 17 weather radars is currently installed, with three more planned to be constructed in the near future (see Fig. 1a).The network was established in 2011 and is managed by the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA), the country's central weather institution.The radar network aims to provide coverage for most parts of Philippines, a country that is frequented by tropical cyclones, monsoon rains, and thunderstorms.One of the network's dual-pol C-band radar is located at the top of the Tagaytay Ridge at the southern portion of Luzon Island.This radar  (hereby referred to as Tagaytay radar) makes a 360-degree scan every 10-min at 16 elevation angles.One commonly used method to retrieve rainfall data from radar reflectivity is the use of the rainfall rate -radar reflectivity [R(Z)] relation (Marshall et al. 1947;Marshall and Palmer 1948;Jorgensen and Willis 1982).The R(Z) relation takes the power form Z = aR b with different regions adapting varying a and b values.The performance of a specific R(Z) relation is influenced by the wide variability in geographic location, season, precipitation phase, intensity, and type of rain event (Rosenfeld et al. 1993;Wu et al. 2018).The application of a single R(Z) relation to regions with different characteristics risks the under-or overestimation of the actual rainfall.One way to calibrate R(Z) relations within a certain region is through the use of raindrop size distribution (DSD) which is the variability of number and diameter of rain drops within a certain volume.Various studies show that DSD measurements around the world vary significantly both temporally and spatially, and can be linked to the type of precipitation, climate, atmospheric conditions, and geographical setting of an area (Lam et al. 2011;Marzuki et al. 2013;Hachani et al. 2017).In retrospect, both Z and R are highly dependent on DSD variability within a radar measurement, as the latter is usually assumed.Marshall and Palmer initially suggested various forms of a and b coefficients and were later complemented and expanded by various authors (Marshall and Palmer 1948;Fulton et al. 1998).Globally, the most commonly used R(Z) relation is Z = 300R 1.4 developed by the United States' National Weather Service (NWS) and for a while was the default R(Z) relations for the United States' WSR-88D radar network (Woodley and Herndon 1970;Kuligowski 1997).As aforementioned, however, R(Z) is best calibrated within a specific region first to improve its performance in rainfall retrieval.
Apart from the use of R(Z) relations, rainfall can also be calculated from dual-pol radars using polarimetric variables such as differential reflectivity (Z DR ), differential phase (Φ DP ), and specific differential phase (K DP ) (Carey et al. 2000;Ryzhkov et al. 2005;Thompson et al. 2018).The use of these dual-pol variables has some advantages over the use of R(Z).For example, K DP is unaffected by partial beam blocking, attenuation, and wet radome effects (Chandrasekar and Cifelli 2012).Moreover, studies found that R(K DP ) usually provide the best rainfall estimation among various QPE relations, as K DP is proportional to the 4.24 th moment of DSD, and is closer to the 3.67 th moment (proportional to rain rate) than the 6 th moment (proportional to Z) (Sachidananda and Zrnić 1986;Ryzhkov et al. 2005).As with the case of R(Z) relations, the use of polarimetric variables can also be improved when their respective relations with R are calibrated on a regional scale (Zhang et al. 2019).Commonly used relations are R(Z, Z DR ), R(K DP ), and R(K DP , Z DR ) expressed as R = aZ b Z DR c , R = aK DP b , and R = aK DP b Z DR c , respectively.Z and Z DR are expressed in mm 6 m -3 , K DP is in ° km -1 , and R is in mm hr -1 while a, b, and c are derived constants (Bringi and Chandrasekar 2001;Montopoli et al. 2017).
In the Philippines, previous studies on polarimetric QPEs use pre-calculated relations derived from other areas (Fulton et al. 1998;Heistermann et al. 2013;Crisologo et al. 2014) and currently, there is a lack of studies aimed at calibrating the dual-pol rainfall relations for any part of the country.This study is hence the first to estimate localized and appropriate R(Z), R(Z, Z dr ), R(K DP ), and R(K DP , Z DR ) relations that best fit the climatic regime of Southern Luzon, Philippines using DSD data.To achieve this objective, the raindrop characteristics of the area of interest was first studied using DSD data from a newly installed optical Parsivel disdrometer located inside the University of the Philippines Los Baños during the wet season of 2018 (see Fig. 1b).The analyzed DSD data was eventually utilized to estimate various dual-pol rainfall relations and was evaluated using four nearby rain gauges within the area with the goal of improving QPEs from the Tagaytay radar.Such DSD analysis, along with DSD-derived radar calibration, is rarely studied in the Philippines.

Study Area, Event of Interest, and Sources of Rainfall Data
The Philippines is characterized by a Tropical monsoon climate according to the Köppen climate classification system.The rainy season from June to October is dominated by the Southwest Monsoon, and is driven by the huge water mass of the Indian Ocean (Flores and Balagot 1969;Cayanan et al. 2011).This period is also where the majority of the country's 18 to 20 tropical cyclones form annually.For this study, the location of the disdrometer, rain gauges, and radar used is around Southern Luzon (Fig. 1b).The study's event of interest is the monsoonal rain event enhanced by Tropical Storm Yagi (hereby termed Yagi from hereon) from 10 to 11 August 2018 UTC.The tropical storm itself formed on 8 August 2018 UTC to the Northeast of Luzon Island with winds reaching up to 35 knots.Although Yagi itself did not make landfall in the Philippines, the storm's peak winds of 75 km hr -1 and continuously intensifying structure enhanced the southwest monsoon that caused massive flooding in Metro Manila and the surrounding provinces.According to the National Disaster Risk Reduction and Management Council (NDRRMC), a total of 187744 families/828462 persons within 680 barangays in five regions (including the National Capital Region) were affected (National Disaster Risk Reduction and Management Council 2018).The monsoonal event brought about by Tropical Storm Yagi was chosen because of the following considerations: (1) concurrent availability of radar data from Tagaytay radar with minimal missing time slots (radar data from October 2018 until the conduct of this study are either missing or unusable due to damaged equipment); (2) availability of rainfall data from four rain gauges near the disdrometer for the whole period; and (3) the event was the longest rainfall event recorded by the Parsivel disdrometer during the rainy season of 2018.Four rain gauges near the disdrometer were chosen as reference data for QPE validation: (1) Agro-meteorological Station where the disdrometer is also located; (2) Institute of Plant Breeding Automatic Weather Station; (3) New College of Arts and Sciences Building Automatic Weather Station; and (4) Makiling Botanic Gardens Automatic Weather Station.The first gauge has a 10-min measurement interval while the other three stations have 15-min time intervals.The observational periods and number of samples from the instruments used in this study is summarized in Table 1.

Raindrop Size Distribution (DSD)
Raindrop size distribution (DSD) is defined as the number of drops within a 3-dimensional space for each diameter class and is usually measured using a disdrometer.Marshall and Palmer initially described DSD as an exponential function given by: where N 0 (the scale parameter) = 8000 and Λ (or the slope parameter) = 4.1R -0.21 (Marshall and Palmer 1948).However, this was developed using data from mid-latitude stratiform rain and was observed to not be able to accurately represent DSDs from other climate types (Williams and Gage 2009).
Ulbrich was the first to propose the use of gamma distribution to provide a better fit especially for very small and very large raindrops (Ulbrich 1983).This form is composed of three parameters: N 0 or the intercept parameter (m -3 mm -1-µ ), n or the shape parameter, and Λ (mm -1 ) or the slope of the distribution and is given by the equation: Other authors such as Zhang et al. (2001) and Brandes et al. (2003) suggested a mathematical relation that exists between n and Λ and further redefined the gamma distribution by assigning specific values to N 0 .Despite the existence of other forms of raindrop size distribution such as the Log-Normal distribution (Feingold and Levin 1986), Scaling Law distribution (Torres et al. 1994), and Weibull distribution (Best 1950), Ulbrich's three-parameter gamma distribution is still the most widely used form by meteorologists and atmospheric physicists for DSD modelling (Ulbrich 1983;Wong and Chidambaram 1985;Smith 2003).
This study utilized the DSD measured by an optical Parsivel disdrometer located at the Agro-meteorological station inside the University of the Philippines Los Baños, Philippines.It uses a laser beam emitter and a mounted receiver, by which the signal voltage measured changes when passed by a hydrometeor.The Parsivel disdrometer measures the particle size and velocity of incoming precipitation simultaneously, groups them into 32 bins (see Table 2), and derive the DSD of a rainfall event using the equation: where N(D i ) is the drop concentration with diameters from D i to D i + ΔD i per unit size interval, n ij is the number of drops within the size bin i and velocity bin j, A (m 2 ) is the sampling area, Δt (s) is the sampling time, D i is the raindrop diameter for the size bin i (in mm), and ΔD i is the corresponding diameter interval (in mm), and V j is the fall speed for the velocity bin j (Löffler-Mang and Joss 2000).Parsivel disdrometers also calculate the type, amount, intensity, and kinetic energy of the precipitation and has the capability to identify the type of precipitation as well as to derive its rainwater content (W), reflectivity (Z), and rainfall intensity (R) from those measurements using the following equations: where W is the rain water content in g•m -3 , R is the rain rate in mm•hr -1 , Z is the radar-reflectivity in mm 6 m -3 , D i 6 is the 6 th power of the diameter of the i-th particle in mm, N(D i ) is the number of particles of diameter D i , and V j is the fall velocity for bin j in m s -1 .In fitting the observed DSD to a gamma distribution model, the parameters N 0 , n, and Λ were first calculated in terms of DSD moments defined by the equation: where M n is the nth moment of N(D) and Γ(x) is the complete gamma function.The three gamma parameters are calculated using the following equations: where: Apart from these integral parameters, the mass-weighted volume diameter (D m ) can also be calculated from DSD using the equation: Moreover, the normalized parameter of N 0 that is not dependent on the parameter n called N w as proposed by Testud et al. (2001) was also calculated using the following equation: where w t is the density of water (1.0 g cm -3 ).Although optical disdrometers have been widely employed by various institutions around the globe because of its portability, quantitatively reliable results, and low cost, there are still some sources of error that may significantly affect its accuracy (Caracciolo et al. 2006;Battaglia et al. 2010;Thurai et al. 2011;Friedrich et al. 2013).For example, Jaffrain and Berne (2011) found out that the total concentration drop uncertainty and radar reflectivity uncertainty can exceed 10% for 1-min samples in Switzerland.Moreover, Krajewski et al. (2006) concluded that the instrument's background noise might be the cause of a higher concentration of smaller drops and higher rainfall rates when compared to other instruments (i.e., 2D video disdrometer and pluviometer).Wind can also heavily affect the reliability of the data measured by an optical disdrometer.Friedrich et al. (2013) found out that unrealistic fall velocities caused by the tilting of raindrops was observed when wind speed exceeds 20 m s -1 and with some occurrence when wind speed is as low as 10 m -1 .Hence, various data quality correction methods and filters must first be employed to make sure that the disdrometer data to be analyzed is sound.
The study utilized the Parsivel DSD measurements from 00:00 UTC 1 June 2018 to 23:59 UTC 31 October 2018 to evaluate the rainfall characteristics of the region during the wet season.Data using 1-, 2-, 5-, and 10-min integration intervals (Δt) were evaluated for the duration of the study.The higher time steps involve the averaging of the 1-min variables over the longer time step (e.g., R).To make sure that the DSD data are sound, the following quality control procedures were applied: (1) Rainfall events characterized with less than 0.1 mm hr -1 and/or having less than 10 droplets were considered noise and were discarded from the analysis; (2) Removal of data when the observed velocity has greater than 40% bias compared to the ideal velocity; and (3) Data are only logged when wind conditions does not exceed 10 m s -1 (Friedrich et al. 2013).
The calculation of radar variables was done using the T-matrix/Mueller method initially proposed by Waterman (1965Waterman ( , 1971)).This method calculates electromagnetic scattering by single non-spherical raindrops and are defined by radio frequency, temperature, type of hydrometeor, canting angle, and axis ratio.The radar specifications used in T-Matrix/Mueller method are shown in Table 3.The assumed axis ratio is from the study by Brandes et al. (2002) from numerical simulations and wind tunnel tests employing a fourth-order polynomial equation: . . .

D D D D
0 9951 0 0251 0 03644 0 005303 0 0002492 where c is the axis ratio and D is the raindrop diameter.
Reflections are assumed to exhibit the "Mirror" reflections symmetry about a plane as defined by Bringi and Chandrasekar (2001).

Radar Data Processing Overview
The study made use of data retrieved from PAGASA's dual-polarized C-Band Tagaytay radar located at the highest peak of the Tagaytay ridge (14.123°N, 120.974°E at 752 masl) in Southern Luzon, Philippines (see Fig. 1).It was installed in April 2012 and originally gathers a full volumetric scan every 15 min and currently scans every 10 min.Table 4 provides the technical specifications of the Tagaytay C-band radar.The radar is about 25 km away from the Parsivel disdrometer, where the lowest available data is at about 1 km above the surface.The dataset used in this study is from 22:00 UTC 10 August 2018 to 16:00 UTC 11 August 2018, where 10-min scans from 09:10 UTC to 10:50 UTC at 11 August 2018 were missing.The processing of radar data is as follows: (1) Clutter detection and removal using co-polar cross correlation coefficient (ρ); (2) Differential phase (Φ DP ) processing and Derivation of Specific Differential Phase (K DP ); (3) Attenuation correction using Specific Differential Phase (K DP ); (4) Systematic bias calculation and correction via weak echo region (for differential reflectivity Z DR ) and using theoretical K DP values (for reflectivity Z HH ); and (5) Interpolation of plan position indicator (PPI) variable values into constant altitude plan position indicator (CAPPI), 3-D gridding, and eventual retrieval of Quantitative Precipitation Estimate (QPE) or rainfall rate (R).Each section's details and specifications are also outlined in this section.

Clutter Detection and Removal Using Co-Polar Cross Correlation Coefficient (ρ)
The removal of non-meteorological echoes is a crucial prerequisite in radar data quality control.In this study, the co-polar cross-correlation coefficient ( HV t or simply t), which is the correlation between the horizontal and vertical co-polar echoes, was used for clutter detection and removal.It is defined as: where V is the complex echo signal sample received from h(v) as the horizontally polarized echo sample, m is the echo sample number, and M is the total number of samples (Li et al. 2014).Figure 2 shows an example image of HV t -based clutter detection and removal dated 08:30 UTC 11 August 2018 using the minimum threshold HV t = 0.95.Notice the removal of non-red signals after filtering.

Differential Phase (Φ DP ) Processing and Derivation of Specific Differential Phase (K DP )
As the differential propagation phase Φ DP tends to increase as the range distance increases, phase folding can occur when it exceeds a certain maximum unambiguous value (in this case is equal to 180°).Hence, unfolding is the first step that was executed in differential phase processing.In this study, the algorithm described by Wang and Chandrasekar (2009) was applied.In addition to unfolding, subtraction of the initial beam phase was also done to remove the system offset.For this study, the initial Φ DP was calculated as the mean Φ DP from each volume scan.Figure 3 shows the results of differential phase processing.Note that the initial values beyond 180° were removed after subtraction.Specific differential phase K DP was derived using the smoothed Φ DP data by least squares fit of multiple differences within 25 bins centered at the range bin.

Attenuation Correction Using Specific Differential
Phase (K DP ) One of the most common methods to calculate polari-metric radar attenuation (A H ) is through the use of the specific differential phase (K DP ) which has a linear proportionality to cumulative attenuation in common radar frequencies (Bringi et al. 1990).Bringi et al. (2001) proposed a selfconsistent method, in which an optimal proportionality of horizontal attenuation is estimated by comparing calculated (or theoretical) K DP with the measured K DP .In this study, the coefficients are derived by fitting lines for A H and A DP (differential attenuation) versus K DP derived from 2018 wet season's DSD data (see Fig. 4).These variables were calculated from DSD using the following equations: where A H,V is the horizontal and vertical attenuation (dB km -1 ), Q H,V (D) are the extinction cross sections for H and V polarized waves, A DP is the differential attenuation (dB km -1 ), K DP is the specific differential phase (in degree km -1 ), λ is the wavelength, and f H and f V are the forward scattering amplitudes for horizontally and vertically polarized waves (Bringi et al. 1990).Figure 4 shows the fitted A H and A DP versus K DP which derived the following relations: . A K 0 0625 Figure 5 shows Z HH images before and after attenuation correction using data from 08:30 UTC 11 August 2018 of the Tagaytay radar.

Systematic Bias Correction
After attenuation correction, the determination and correction of systematic bias was executed for differential reflectivity (Z DR ) and reflectivity (Z HH ).Z DR bias is usually determined through the use of a vertically pointing radar, wherein the shape of the raindrops as seen at 90° elevation angle is assumed to be nearly circular (Gorgucci et al. 1995;Vivekanandan et al. 2003).In the absence of vertical measurements as is the case of this study, the use of reflectivity data from light rain is suggested by Smyth and Illingworth (1998) who pointed out that rain drops in very light rainfall regions are almost spherical.In order to calculate the bias, gates with the following characteristics were utilized and classified as light rain:   (1) ΔΦ DP < 10° km -1 ; (2) 10 dBZ < Z HH < 20 dBZ; (3) ρ HV > 0.98; and (4) altitude < 3.5 km (to avoid freezing height).
The value of Z DR bias was calculated as the average of derived bias values of all beams from all elevation angles from each radar sample as shown in supplementary Fig. 1a.The calculated mean Z DR bias value for all samples is 0.61 dB.
For Z HH bias calculation, the self-consistency approach proposed by Scarchilli et al. (1996) wherein Z HH bias is defined as the discrepancy between calculated and theoretical K DP values and is calculated as: For this study, DSD-derived Z HH (mm 6 m -3 ) and Z DR (linear units) data were fitted with K DP (deg km -1 ) for the whole 2018 wet season which resulted to the following model calculations: c m in dBZ.Z HH bias was calculated for each radar sample as difference between the first 5 and last 5 gates and were averaged over all elevation angles and azimuthal angles.Supplementary Fig. S1b shows the different Z HH bias values for each radar sample, as well as the final Z HH bias (-0.333 dB) which is calculated as the mean of all the Z HH bias values.
Both the calculated Z DR and Z HH bias values (0.61 and -0.333 dB, respectively) are significantly smaller than the values obtained by one of the recent studies of Crisologo et al. (2014) using 2012 data from the Tagaytay C-Band Radar (Z DR bias = 1.2 dB and Z HH bias = -8.5 dB).However, it must be noted that the radar has undergone several damage, repair, and calibration cycles over the years which may have positively affected the quality of the radar data itself.

Gridding and Parameter Retrieval
To obtain gridded values, data were interpolated from Plan Position Indicator (PPI) spherical coordinates into constant altitude plan position indicator (CAPPI) 3-dimensional Cartesian coordinates with a horizontal resolution of 1 × 1 km and a vertical resolution of 0.5 km.Rainfall rate (R) in mm hr -1 was obtained using the power law Z = aR b where a and b were set from either: (1) standard values from literature; or (2) obtained from this study (see section 2.5).All R values were calculated in spherical form before gridding.The final R value for a certain point was calculated as the median value between a 3 × 3 grid-area around the location of interest.

Estimation of R(Z), R-Z-Z DR , R-K DP , and R-K DP -Z DR
Relation Parameters R(Z, Z DR ), R(K DP ), and R(K DP , Z DR ) relations were calculated using least-squares nonlinear fitting of 1-min data from the Parsivel disdrometer.On the other hand, two methods were used to calculate R(Z) relations in the form Z = aR b .The first method (herein named as Method 1) is the use of DSD gamma distribution while the second method (herein named as Method 2) utilizes linear fitting of log-transformed R(Z) pairs to extract a and b values.For Method 1, The second method includes fitting a linear regression equation between R(Z) pairs.To fit their power law form, however, the values were first converted to their natural logarithmic forms.From the equation: where R is in in mm hr -1 and Z is in mm 6 m -3 , Equation ( 29) can be further simplified by assigning new variables: . Z is then expressed to Z dBZ using the following: From Eq. ( 34), the following equation to calculate R was generated: which is linear in nature.Hence, simple linear regression (Freedman 2009) was applied and a and b were extracted from the resulting linear equation using the following equations: a e e As with the first method, the natural logarithms of a values were also first retrieved to filter unrealistic values.

Evaluation of the Derived Relations
Rainfall values derived from various relations and aggregated into 30-min interval time step were compared using six statistical validation variables: Pearson's correlation coefficient (r), mean error (ME), Nash-Sutcliffe efficiency (NSE), mean absolute error (MAE), percent bias (pBIAS), and root-mean-square error (RMSE).Pearson's r measures the goodness of fit and linear association between the QPE (P) and the gauge data (G), and is given by the following equation: where P and G are the mean values of the QPE product and gauge data.Values of r range from -1 to 1, where a positive r value denotes positive correlation while a negative r value implies negative correlation between the datasets (Pearson 1895).A zero r denotes that the datasets exhibit no linear correlation while a -1 or 1 r value means a perfect linear relationship.Mean error (ME) is the mean value of the difference between the QPE (P) and the gauge data (G), as given by the following equation: A zero ME value denotes the absence of error.A negative value implies underestimation while a positive value means an overestimation of the observed dataset.The absolute value of ME is given by the mean absolute error (MAE): Percent bias (pBIAS) measures the mean error of the observed data against the reference data expressed in percentage and is computed as: The Nash-Sutcliffe efficiency (NSE) is a normalized statistic that determines the magnitude of the residual variance in comparison with measured data variance (Nash and Sutcliffe 1970).NSE is calculated using the following equation: where NSE can be anywhere within the range (-inf, 1).An NSE value of zero indicates a perfect match of the datasets without any errors while a value of 1 indicates that the modelled data is as accurate as the mean of the observed data.
Values between 0 and 1 denote some deviations from the observed dataset.A negative NSE value indicates poor quality of the modeled data and implies that it is not fit as an estimator for the observed dataset (Grunwald and Frede 1999).
The root-mean-square error (RMSE) is the square root of the mean squared differences (errors) between QPE (P) and gauge data (G), and is given by the following equation:

DSD-Derived Rainfall (R DSD )
Figure 6a shows the relationship between observed rainfall from rain gauge (R RG ), and DSD-derived rainfall (R DSD ) in mm.A total of 1126 pairs of co-available 10-min datapoints were used for this plot.On the other hand, Fig. 6b shows the comparison between the accumulated rain for the whole wet season (mm), respectively.The observed total rainfall for the whole period is 833 and 816.12 mm for rain gauge and Parsivel, respectively.These values agree with the report of Bagtasa (2019), wherein the authors found out that out of the approximately 2000 mm of annual rain fall within the Philippines, 43% (860 mm) is derived from the southwest monsoon.Results show that a generally good agreement is exhibited by R DSD and R RG (r = 0.86).However, this correlation is comparatively low and can be attributed to the low data resolution of rain gauge (with a minimum detection of 0.5 mm in comparison to modern rain gauges with 0.1 to 0.25 mm limit of detection).In addition, there are also some datapoints where the rain gauge registered zero rainfall while the disdrometer measured high rainfall, which can be caused by instrument error in the tipping bucket rain gauge.Nevertheless, the accumulated rainfall data shows very high agreement, with only 2% (16.88 mm) discrepancy and suggests that the data derived from the Parsivel disdrometer can be utilized for further processing.The composition of the DSD-derived rainfall products was also analyzed for the whole period and for Yagi (supplementary Fig. S2 shows exemplary R distribution in terms of 1-min Δt) while Table 5 shows the general statistics of the datasets at various integration time steps.Results show that the wet season has an average of 4.08 to 4.79 mm hr -1 rainfall rate for different time integration steps.From the table, it can also be observed that Yagi has less mean and maximum R values than the whole of wet season, which states that the general rainfall rate in this particular event is less than that of the season it belongs to.Overall, 82% of all the DSDderived rainfall data throughout the wet season is characterized by rainfall events below 5 mm hr -1 for 1-, 2-, and 5-min Δt and 83% for 10-min Δt.

Drop Size Distribution and Normalized Gamma Parameters
Figure 7 shows the average raindrop diameter (D in mm) for the whole wet season and for Yagi at various integration time steps against the average drop number [N(D) in m -3 mm -1 ] in log-scale.Results show that the data for wet season at 1-min Δt exhibits the widest range of observed raindrop diameter with samples exhibiting diameters beyond 8 mm, as well as the highest average number of drops for each bin size.On the other hand, the samples observed for Yagi for all integration time steps exhibit smaller average diameters and fewer number of drops than the overall mean.It can be observed from the plot that the number of raindrops generally increases and have peaks at around 1 mm diameter and then gradually decreases from thereon.These are the same results as with the study in some neighboring regions such as Singapore (Kumar et al. 2010), Malaysia (Lam et al. 2015), and East China (Chen et al. 2013).It can also be observed in the plot that the short time steps include larger rain drops than the long time steps, which can be attributed to the averaging effect of using longer measuring periods.
The general statistics of the calculated normalized gamma parameters (D m , log 10 N w , n, and Λ) are summarized in Table 6.Skewness and kurtosis are calculated respectively as:   and Yagi (broken lines) using 1-, 2-, 5-, and 10-min Δt.

SK E x x
where SK is skewness, x is the sample, x is the sample mean, E is the expectation of (xμ), KT is kurtosis, and N is the sample number.In addition, the histograms of D m (blue, in mm) and log 10 N w (red, in mm -1 m -3 ) are shown for the wet season and Yagi using all four integration time steps in Figs.8a and b.Lastly, the values of the gamma parameters were compared with the results of studies conducted by some neighboring countries (see Table 7).In terms of D m , results show that the average (1.36 to 1.48 mm) and standard deviation (0.31 to 0.46 mm) values are close to each other for the whole wet season and Yagi on all time integration values.Additionally, all datasets have positive skewness (or skewed to the left) which denotes that the number of occurring small raindrops are higher than the large ones.Positive kurtosis values for all datasets show that the number of occurrences of raindrop D m values are near the overall mean.Comparing the mean D m with other studies, results show that the mean raindrop diameter in the   Bringi et al. (2003Bringi et al. ( , 2009)).Note: 1: not stated in their respective published studies.

Country
Philippines for the wet season (at 1-min Δt) is the same as the value observed by Tokay and Short (1996) for the tropical ocean (1.41 mm).However, the mean D m in this study is observed to be lower than the study by Lam et al. (2015) in Malaysia (1.74 mm), and from the results of Chen et al. (2013) in Eastern China (1.66 mm), but is observably higher than that from the study of Lee et al. (2019b) in Northern Taiwan (1.16 mm), from the study by Seela et al. (2017) in Palau (1.11 mm) and from 1-min DSD results in Taiwan from the same study (1.24 mm).These results imply that raindrops in the Southern Luzon, Philippines is generally smaller than that from Malaysia and Eastern China but are larger than those from Taiwan and Palau.For log 10 N w values, the mean and standard deviation values from both periods of all integration time steps have minimal differences (2.94 to 3.18 m -3 mm -1 for mean, 0.45 to 0.55 m -3 mm -1 for SD).Significant differences, however, was observed for the two dataset periods in terms of skewness.For the whole wet season, the log 10 N w values tend to be negatively skewed (the frequency of occurrence is toward the higher values) while for Yagi, they were observed to be positively skewed (the frequency of occurrence is toward the higher values).This denotes that there are fewer number of raindrops for each class in Yagi than for the wet season in general.It was noted that the maximum log 10 N w values for the whole wet season (3.18 m -3 mm -1 at 1-min Δt) is observably less than that from Taiwan (4.22 and 4.36 m -3 mm -1 for 1-and 10min, respectively) as well as Palau (4.56 m -3 mm -1 ), and less by a small amount in comparison with that from Malaysia (3.52 m -3 mm -1 ), and East China (3.42 m -3 mm -1 ).This implies that compared to these areas, the number of raindrops is less in Southern Luzon, Philippines than its neighboring regions.It must be noted, however, that the mean R from this study is observably less than the aforementioned studies (see Table 7).The shape parameter µ gives an idea of how spread the distribution of raindrops is across the DSD classes.A big difference was observed between the µ values from the 1-, 2-, 5-, and 10-min integration time steps for all time periods.The value of µ is significantly higher for 1-min integration time step (11.20 for wet season, 9.52 for Yagi) than for 10min integration time step data (6.47 for wet season, 5.39 for Yagi), with in-between values using 2-and 5-min integration time steps.This denotes that the data from 1-min integration time steps covers a wider range of size distribution than the data from 10-min integration time step.However, a reason for such high distribution might be the presence of outlying number of drops on larger diameters that affects the value of µ.Moreover, the average µ value for 1-min data for the wet season (µ = 11.20) is also significantly higher than the µ values of China (µ = 3.51, 1-min Δt), Palau (µ = 8.37, 1-min Δt), Taiwan (µ = 6.72 and 5.60 for 1-and 10-min Δt, respectively), and Malaysia (µ = 6.76, 1-min Δt) (Chen et al. 2013;Lam et al. 2015;Seela et al. 2017;Lee et al. 2019b).However, when using longer time-steps (> 1-min), the values of µ appears to be closer to the aforementioned areas of comparison, albeit still a little higher in magnitude.Nevertheless, it can be concluded that rainfall in South Luzon is more widely distributed across DSD classes than the aforementioned neighboring areas.Significant differences were also observed between the four integration time steps in terms of the slope parameter Λ (mm -1 ).In gamma distribution, a higher slope parameter signifies a gentler slope of the gamma curve.Data using 1-min Δt (wet season = 13.04 mm -1 , Yagi = 11.05mm -1 ) achieved significantly higher Λ values compared to those using 10-min Δt (wet season = 8.23 mm -1 , Yagi = 7.05 mm -1 ), with in-between values using 2-and 5-min integration time steps.Compared to nearby areas, the achieved Λ using 1-min Δt from the study area is also significantly higher than those from China (Λ = 4.52 mm -1 ), Malaysia (Λ = 7.34 mm -1 ), and Taiwan (Λ = 9.82 mm -1 ) but is close to the results in Palau (Λ = 12.19 mm -1 ).As with the case of µ, the use of longer time steps achieved Λ values closer to that of its neighbors.
The scatter plots of the normalized parameters log 10 N w (in mm -1 m -3 ) and D m (in mm) are shown in Figs.8c and d.Comparing the scatter plots of wet season and Yagi for both integration time reveals a generally common mean, with the major difference only on the number of samples plotted.The blue and red squares in the plot are the range of maritime convection and continental convection (respectively) as reported by Bringi et al. (2003Bringi et al. ( , 2009) ) by studying DSDs across various climatic regimes.In addition, the study also reported the separation between convective and stratiform rain DSDs which is shown in dashed line in the plots.Results of the study showed that the mean and most of the plotted data are overwhelmingly on the stratiform-class side of the dashed line, which denotes that most of the rainfall that occurred during the wet season and Yagi are stratiform in nature.It was also observed that in Southern Luzon, Philippines the convective rains are closer to continental than maritime in terms of origin.
To further investigate the variability of log10(N w ) and D m with respect to rain rates, scatter plots of the two parameters vs R (mm hr -1 ) are also shown in Fig. 9.In addition, the plots also show the fitted power law curve and relationship equation derived using least-squares method.A general positive trend of the curve's exponent can be observed for all the plotted data in terms of log10(N w ) and D m .This denotes that for all the datasets, both parameters increase with heavier rainfall intensity, owing to more efficient coalescence and breakup mechanisms (Hu and Srivastava 1995).In terms of log10(N w ), the coefficients of datasets using 1-min time integration steps (3.1369 for wet season, 3.0663 for Yagi) are higher than the datasets utilizing 10min integration time steps (3.0112 for wet season, 2.9792 for Yagi), with in-between magnitudes using 2-and 5-min integration time steps.However, the inverse is observed in the curve exponents, wherein 1-min integration time steps have lower exponent values (0.054 for wet season, 0.071 for Yagi) than 10-min integration time steps exponents (0.0627 for wet season, 0.0774 for Yagi).This denotes that although log10(N w ) in 1-min Δt have a higher intercept, the rate of growth as R increases is lower compared to that of 10-min Δt.In the case of R vs D m , however, all coefficients and their respective exponents are relatively close to each other in terms of integration time steps.

DSD-Derived Reflectivity (Z DSD )
The histograms of the calculated reflectivity Z values in dBZ for are shown in Fig. 10.Results show that for the wet season, most of the measured reflectivities are mostly within the 20 -30 dBZ range for the data using 1-, 2-, and 5-min Δt (figures are not shown for 2-and 5-min data).On the other hand, data using 10-min Δt has observably lower measured Z, with most of the data included in the 18 -22 dBZ range (figure for 10-min is not shown).In terms of the results for Yagi, the histograms for all integration time steps show that most of the data are within the 18 -25 dBZ range, and are hence considered to have lower reflectivities than the overall average for the wet season.It was also observed that the data have close kurtosis values (with a range of 3.04 to 3.35) and are all observably skewed to the left (positive skewness), which indicates that the majority of the data are within the lower range of Z.To illustrate the relationship between R DSD and Z DSD , scatter plots for log 10 R (mm hr -1 ) vs. Z (dBZ) are shown in Fig. 11.Included in the plots are the mean (black dot) and the standard deviations (black lines propagating from the black dot) of the data.Also shown in the plots is the separation line between convective and stratiform rain events in terms of R and Z by Bringi et al. (2003Bringi et al. ( , 2009)).Results show that an overwhelming majority of the rain events measured by the Parsivel were of stratiform origin in terms of both R and Z.This pattern can be observed in all datasets during the wet season and for Yagi, and for both 1-and 10-min integration time steps.These results agree with the result observed from the D m vs. log 10 N w plots shown in Figs.8c and d.

Radar Reflectivity Retrieval
The time-series comparison of Z DSD and reflectivity measured by the Tagaytay radar (Z RADAR ) is shown in Fig. 12a while the scatter plot is shown in Fig. 12b.Results show varying trends in the agreement between Z RADAR and Z DSD with respect to its magnitude.For example, between 22:00 UTC 10 August and 00:00 UTC 11 August, as well as between 11:00 and 16:00 UTC 11 August, the agreement between the two datasets is generally good, with values ranging from 15 to 33 dBZ.However, the discrepancy between the two datasets increases towards the main rainfall event, which is the period between 4:00 and 10:00 UTC 11 August (where the peak rainfall for Yagi occurs).Although the patterns of peaks and troughs of the time series plots for both datasets are still comparable, their differences are more pronounced in this period, with the highest radar bias of 15.87 dBZ at 5:30 UTC 11 August.Nevertheless, in the period with the highest burst rainfall (between 08:00 and 09:00 UTC 11 August), the two datasets have a generally good agreement and minimal bias values.Overall, the Z DSD provided an acceptable agreement with the Z RADAR and were hence used to calculate QPEs using different R(Z) relations.

Derived Polarimetric Rainfall Relations
The generated a and b values for R(Z) relations are summarized in Table 8.The R(Z) Code column refers to the  codename given to each R(Z) relation and shall be used to refer to the 16 generated relations from hereon.In general, the values of a using Method 1 (339 to 421 for wet season, 363 to 453 for Yagi) are higher than those from Method 2 (251 to 332 for wet season, 297 to 355 for Yagi) for both time periods.However, the values of b appear to have an opposite pattern as results show that b values using Method 1 (1.19 to 1.24 for wet season, 1.19 to 1.24 for Yagi) are smaller compared to b values derived from Method 2 (1.5 to 1.54 for wet season, 1.34 to 1.41 for Yagi).The final values of a ranges from 252 (A5) to 421 (A4) for the wet season while for Yagi, the final values are between 297 (Y5) and 453 (Y4).It can be observed that Method 1 produces the maximum a values for both time periods, which indicate that Method 1 is more sensitive to the presence of larger rain drops in its R(Z) relations compared to that of Method 2. For parameter b, the final values are between 1.19 (A1) and 1.54 (A5) for wet season, and between 1.19 (Y1) and 1.41 (Y5) for Yagi.
The calculated values of a and b from the gamma parameters µ and N 0 are shown in Fig. 13.The black broken lines represent the mean value of the plot that are used as final coefficients and exponents.For the whole wet season, results show a decreasing variability in magnitudes of a as the time increment increases while values of b have no discerning pattern.This also holds true for the case of Yagi.In terms of magnitude, the values of a and b increases as time increment increases for both the wet season and for Yagi.The occurrence of a constant pattern for magnitudes of both a and b can be attributed to the fact that both variables are calculated as a function of µ, which has shown to have a consistent decreasing pattern as the time integration step increases.For linear fitting, the scatter plots of ln10/10 × Z (dBZ) and lnR (mm hr -1 ) are shown in Fig. 14 for both time periods using all four integration time steps.Results show similar trends with the a values from Method 1, wherein the magnitude of a values increases as Δt increases for both the wet season and Yagi.However, a main difference from Method 1 is that there is an inverse relationship for the observed magnitudes of b and time increments when using Method 2 for both wet season and Yagi.One important observation from the derived R(Z) linear fits is that the residual variance remains stable as time increment increases.The same results were observed by Chapon et al. (2008) from their study of Mediterranean rainfall in France.According to the authors, such a pattern is unexpected given that time integration from 1-to 10-min has the potential to bring sample noise reduction from the data.
The generated R(Z, Z DR ), R(K DP ), and R(K DP , Z DR ) relations from wet season and for Yagi (accompanied by a subscript A and Y, respectively) are also shown in Table 8, while the fitted lines are shown in Fig. 15.Results show that the coefficient a of R(K DP ) and R(K DP , Z DR ) tend to be larger for the Yagi event than the wet season while it is the op-posite in terms of R(Z, Z DR ).As the DSD values were used for the derivation of these relations, the results indicate that Yagi has less concentration of small raindrops and a lower proportion of large raindrops because the coefficient values of a is larger for R(K DP ) and R(K DP , Z DR ).On the other hand, R(Z, Z DR ) relations show Yagi has a lower proportion of large raindrops but a higher concentration of small raindrops, which is why the coefficient a is smaller than the in the wet season.The results implies that rainfall based on R(Z, Z DR ) is more sensitive to the number of small raindrops than the proportion of large raindrops.

Initial Comparison Between R DSD , R RADAR , and R RG during Yagi
To determine the variability between the datasets derived from the three instruments used in this study (rain gauge, disdrometer, and radar), their respective datasets at 10-min time interval during Yagi were first compared.Figure 16 shows the time series plots of R derived from rain gauge, disdrometer, and from Tagaytay Radar using the standard WSR-88D relation Z = 300R 1.4 relation (will be referred to as SR from hereon).In addition, the standard R(Z) relation for neighboring Taiwan from the study of Seela et al. (2017) (Z = 283.35R 1.35 and will be referred to as ST) was also compared to the results of the study.Results show generally similar peaks and troughs between the datasets throughout the study period.However, the discrepancies between the three instruments are more pronounced during peak periods (between 05:00 and 09:00 UTC 11 August) with a maximum absolute bias of 11.88 mm hr -1 (+) for R DSD (at 08:30 UTC 11 August), 13.68 mm hr -1 (-) for R RADAR-SR (at 08:10 UTC 11 August), and 13.02 mm hr -1 for R RADAR-ST (at 08:10 UTC 11 August).In terms of linear regression, the fitted lines show a higher average R for disdrometer compared to that of rain gauge (r = 0.97).On the other hand, an even higher deviation from the rain gauge data was observed for the R RADAR-SR and R RADAR-ST (both have an approximate r of 0.85), albeit unlike the R DSD , both radar products are observed to be characterized by underestimation.The underestimation of R RADAR and overestimation of R DSD is further emphasized through the time series of accumulated rain plots.In comparison with the gaugemeasured total Yagi rainfall (23 mm), R DSD overestimated it by 39% (32 or 9 mm difference) while the R RADAR-SR underestimated the total rainfall by 26% (17 mm) and R RADAR-ST by 27% (18 mm).
The overestimation of DSD rain rate from an optical Parsivel disdrometer has also been observed in other studies (Lanza and Vuerich 2009;Thurai et al. 2011;Tokay et al. 2013Tokay et al. , 2014;;Zhang et al. 2015).The authors of these various studies all found out that the optical Parsivel  disdrometer has a tendency to underestimate smaller rain drops and overestimate larger ones, even in solid forms of precipitation (Battaglia et al. 2010).This also explains why the R derived from the Parsivel disdrometer exhibited a significant amount of overestimation for Yagi (with multiple rainfall bursts) but only a minimal overestimation for the whole duration of the study (which is characterized by moderate to high rainfall stratiform rainfall).Moreover, it was also observed that there is a wide range of raindrop size diameters that is characteristic of the study area as compared to other nearby countries such as Malaysia (Lam et al. 2015) and Taiwan (Lee et al. 2019b).This presence of larger raindrops can further exacerbate the overestimation of the Parsivel disdrometer.In the case of Radar QPE, using the default WSR-88D relation Z = 300R 1.4 was shown to generally create an underestimation within areas outside the area of origin Miami, Florida (Vieux and Bedient 1998;Ulbrich and Lee 1999;Jayakrishnan et al. 2004;Yang et al. 2016).Further analysis by Legates (2000) suggests that this standard R(Z) relation tend to overestimate light rainfall and underestimate heavy rainfall, which can explain why Yagi (which is characterized by multiple bursts of heavy rainfall) was significantly underestimated in this study.As R(Z) relations can vary significantly over a number of factors such as prevailing climate, geography, and other characteristics of a certain region, this trend of underestimation through the use of the standard R(Z) relation is proof that its use is not suitable for Southern Luzon, Philippines.

Statistical Validation of Generated QPEs Using Various R(Z) Relations
As the robustness of a specific R(Z) relation cannot be defined by linear fit and difference in accumulated rainfall alone, six statistical validation tools (Pearson's correlation coefficient or r, mean error or ME, percent bias or pBIAS, Nash-Sutcliffe efficiency coefficient or NSE, mean absolute error or MAE, and root-mean-square error or RMSE) were also calculated for each derived QPE dataset to analyze their respective performances with reference to rain gauge data.Exact values of all the statistical variables from all QPEs and stations are shown in Table 9.There is a highly variable set of results of QPE statistical validation in terms of the six statistical tests applied.QPEs calculated from Method 2 outperformed those from Method 1 in terms of linear relationship r.Chapon et al. (2008) also reported the same results in his study, and stated that it was only natural as regression techniques work directly with DSD moments of interest For r, there is no observed trend when it comes to time integration step and as expected, the highest linear correlation values were observed in the collocated rainfall station, theoretically because both Z and R were calibrated in that area.
For ME, all QPEs were observed to underestimate the actual rainfall and hence were all negative values.An interesting observation is that the two parameters increase in magnitude as the time integration step increases, a pattern which holds true for both methods and from all stations.This can be explained by the increasing magnitude of accumulated rain if the time step is longer.In terms of ME and pBIAS, ST performed better than SR by a significant amount for each individual station and all stations combined.Additionally, Method 1 overwhelmingly outperforms Method 2 for all individual stations as well as the general analysis of all four stations, and 5 QPEs derived from Method 1 were able to outperform SR.In addition to ME and pBIAS, MAE is also a measure of the actual error and is closely related to the two aforementioned parameters.Hence, the results were highly similar with that from ME and pBIAS.MAE was also observed to increase as time integration step increases -potentially because of the same reasons as ME and MAE.
In addition, Method 1 also outperformed Method 2 for all individual stations and as a whole.Five QPEs from Method 1 also were able to outperform SR in terms of ME and MAE.Results of NSE analysis show highly similar trends as the ME, MAE, and PBIAS, as expected because all these variables measure the degree of deviation/error from the observed data.Results show that ST achieved the best NSE value (0.678), followed by A3 (0.673) and Y2 (0.671).This is also the same result that was observed in terms of RMSE, and denotes that ST has the capability to perform well in radar QPE in terms of NSE and RMSE.In comparison to the study by Crisologo et al. (2014) for the Tagaytay Radar using the R(Z) relation Z = 250R 1.25 for hourly rain accumulation for the wet season of 2012, the top-performing QPEs of this study yielded significantly better results in terms of all six statistics and suggests that the derived R(Z) relations has the potential to outperform commonly used relations for radar QPEs (see Table 9, last row).To summarize the results, the QPEs from four derived R(Z) relations which performed consistently well in terms of statistical validation are: A1 (Z = 339R 1.19 ), A2 (Z = 361R 1.20 ), Y1 (Z = 363R 1.19 ), and Y2 (Z = 384R 1.20 ).It must be noted however, that the best QPEs in terms of r was not considered because of very low discrepancies between the r values for all QPEs.SR was not able to perform well at all in terms of any statistics for all the station combined.In fact, even ST outperformed SR for all but 1 statistic (r).This denotes that SR wasn't able to acceptably model the real rainfall for this particular event.

Comparison Between Generated QPEs using Dual-Pol Relations
The best performing R(Z) relations (A1, A2, Y1, and Y2) were compared with the other derived dual-pol relations to determine the most useful dual-pol relation/s for the study area, and the results are shown in Table 10.For this section of the discussion, the aforementioned relations will be referred to as R(Z) A1 , R(Z) A2 , R(Z) Y1 , and R(Z) Y2 , respectively.Additionally, the dual-pol relations have the subscripts A and Y which corresponds to wet season and Yagi, respectively.As shown in the previous section, R(Z) relations underestimated rainfall during Yagi, with varying degrees depending on the specific relation used.When compared with other dual-pol relations, results show that the use of R(K DP ) further reduced this underestimation.In fact, R(K DP ) values outperforms all the other QPEs by a high degree in terms of ME [-0.106 mm hr -1 for R(K DP ) A and -0.133 mm hr -1 for R(K DP ) Y ], MAE [-0.106 mm hr -1 for R(K DP ) A and 0.133 mm hr -1 for R(K DP ) Y ], and pBIAS [-8.314% for R(K DP ) A and -10.400 mm hr -1 for R(K DP ) Y ].The results in the study agrees with various dual-pol QPE comparative analyses in literature (e.g., Crisologo et al. 2014;Chen et al. 2017;Gou et al. 2018;Voormansik et al. 2020) and can be attributed to K DP 's immunity to radar calibration and attenuation (Chandrasekar and Cifelli 2012).However, R(Z) relations, were still able to outperform the others in terms of NSE [0.665 for R(Z) Y1 and 0.658 for R(Z) Y2 ], RMSE [1.809 for R(Z) A1 and 1.781 for R(Z) A2 ], and r [0.842 for R(Z) A2 and 0.839 for R(Z) Y2 ], although the difference is only minimal compared with R(K DP ) in terms of NSE, and with R(Z, Z DR ) in terms of r.Other studies such as Gou et al. (2018) and Voormansik et al. (2020) show that R(Z) has the potential to outperform other dual-pol relations with proper bias and attenuation correction.On the other hand, R(Z, Z DR ) performed poorly for ME, pBIAS, and MAE while R(K DP , Z DR ) performed poorly in terms of NSE and RMSE.This can be attributed to the residual effects such as resonance and noise to Z DR values, which explains why the combination of Z DR with Z and K DP resulted to lower accuracy (Ryzhkov and Zrnić 2019).As the equation constants derived from wet season data and Yagi data have close values, results show that there are minimal differences with their R values (< 5% difference in statistics).Major differences (20 -21%) were observed only in terms of ME, MAE and pBIAS for R(Z) A1, Y1 , and R(K DP ) A, Y , wherein the use of wet season-based R performed better than Yagi-based R. Overall, results show that R(K DP ) performed the best in representing the rainfall of Yagi event when compared to R(Z), R(Z, Z DR ), and R(K DP , Z DR ).

DISCUSSION AND CONCLUSIONS
The present study aims to derive robust and suitable dual-pol relations that can provide acceptable agreement with reference gauge data for the wet season in the Philippines, a country frequented by tropical storms and flooding.The study utilized the DSD data measured by the Parsivel disdrometer at Southern Luzon, Philippines to derive paired rainfall rate -radar reflectivity data using the T-matrix/ Mueller method, and reflectivity data from the Tagaytay Radar to analyze the Yagi Event Monsoon (Yagi) in terms of rainfall and DSD at ground level.The effects of time integration step in DSD data were also analyzed by using 1-, 2-, 5-, and 10-min Δt as inputs for deriving R(Z) relations.Results show that the area's DSD demonstrates a highly variable distribution in rain drop diameter size and is further confirmed by the relatively higher shape parameter µ in its gamma distribution when compared to neighboring countries such as Malaysia, China, Taiwan and Palau (Chen et al. 2013;Lam et al. 2015;Seela et al. 2017;Lee et al. 2019b).It was also observed that on the basis of the normalized DSD parameters N w and D m , the area of interest has a relatively larger average raindrop diameters than its neighboring Asian counterparts, albeit a smaller number in the total number of raindrops in comparison with the same areas.This is supported by the derived slope parameter Λ which is observed to exhibit higher magnitude (flatter slope which denotes a wider truncated gamma curve) than other neighboring countries.R(Z) parameters a and b were derived using three sets of partitions: covered period (wet season and Yagi); time integration step (1-, 2-, 5-, and 10-min); and method used (Method 1 and Method 2) which generated a total of 16 R(Z) relations.Method 1 makes use of the relationship between the gamma parameters µ and N 0 to calculate a and b for each time step while Method 2 derives the two parameters by linear fitting the log-transformed Z and R datasets.Results show that there is a high variability between the calculated a and b in terms of all three partitions.For method type, Method 1 derived relatively higher a values than Method 2, which is likely due to the sensitivity of Method 1 to larger rain drops (as a is calculated from each and every time step).On the other hand, Method 2 utilizes the least squares method, which reduces its sensitivity to outliers in drop size data.However, the calculated b values exhibit an opposite trend, wherein Method 2 consistently achieved higher b values for all combinations that those from Method 1.There is also an observable pattern in a and b values in terms of Δt used.
Values of coefficient a increase as the time integration step increases, which can be due to the inclusion of larger rain drops in a greater proportion of DSD data in longer time steps than in shorter ones.On the other hand, b increases as time step increases for Method 1, while the opposite is observed for Method 2. Lastly, the a values derived from the wet season were relatively smaller than the ones derived from Yagi, which can be attributed to the fact that all weak and strong rainfall events were considered during the wet season while Yagi is considered as a relatively strong event and is therefore more inclusive of larger raindrops.An interesting result on this section is that there is no increase in linear agreement between Z and R as Δt increases, which implies that the aggregation of datasets into larger values does not improve the quality of data despite its expected noise reduction effect (Chapon et al. 2008).
Six validation statistics were used to analyze the performance of the 16 derived QPEs: Pearson's correlation coefficient (r); mean error (ME); percent bias (pBIAS); Nash-Sutcliffe efficiency (NSE); mean absolute error (MAE); and root-mean-square error (RMSE).Besides the 16 derived QPEs, the standard relation Z = 300R 1.4 (SR) and Taiwan's relation Z = 283.35R 1.35 (ST) were also used and compared with reference rain gauge (RG) data.Results for QPE statistical validation varies greatly in terms of the Method used, time integration step, temporal period, and the statistic itself.In terms of simple error statistics (ME, MAE, and pBIAS), Method 1 overwhelmingly outperformed Method 2 for individual stations and for all stations analyzed as a whole.In fact, all five QPEs that were able to outperform SR were derived using Method 1 (with A1 and A2 consistently performed best in terms of these three statistics).However, these errors appear to grow larger in magnitude as Δt increases which is likely due to the higher aggregated rainfall caused by longer aggregation time.On the other hand, Method 2 performed better in terms of linear agreement r, which implies that the general trend of Method 2 mirrors the linear trends of RG better than Method 1. Being based on linear regression itself, this performance of Method 2 can be attributed to the fact that regression techniques work directly with the DSD moments of interest (Chapon et al. 2008).It is important to note that SR was not able to perform well in terms of all computed statistics, while ST was able to rank first in terms of NSE and RMSE.Four R(Z) relations consistently performed well for all the statistics: A1 (Z = 339R 1.19 ), A2 (Z = 361R 1.20 ), Y1 (Z = 363R 1.19 ), and Y2 (Z = 384R 1.20 ).These relations were able to reduce the QPE bias by up to 61% when compared to SR.In addition to the determination of best R(Z) relations, the study also generated R(K DP ), R(Z, Z DR ), and R(K DP , Z DR ) relations and analyzed their intercomparison.Results show that R(K DP ) produced the best results in terms of ME, MAE, and pBIAS.The use of R(K DP ) was also shown to exhibit the most accuracy in many studies (e.g., Crisologo et al. 2014;Chen et al. 2017;Gou et al. 2018;Voormansik et al. 2020) and is usually attributed to the immunity of K DP to radar calibration and attenuation.Nevertheless, the R(Z) relations were still able to perform well in terms of r, NSE, and RMSE.As studies such as Gou et al. (2018) and Voormansik et al. (2020) stated, R(Z) has the potential to outperform other dual-pol relations when bias and attenuation are properly corrected.
The results of the study's intercomparison between standard relations and the derived relations present a significant improvement over the Philippines' general rainfall weather radar rainfall retrieval from single-and dual-pol radars alike.This significance is further emphasized by the strategic location of Tagaytay radar, which is according to Crisologo et al. (2014): (1) provides coverage for the Bataan Peninsula, which is covered by an S-band radar that is severely affected by beam blockage; (2) covers the area of Metropolitan Manila's basins, the country's densest region; and (3) is the only source of data for a number of provinces in Southern Luzon.Additionally, the comprehensive description of the region's DSD also furthers the understanding of the country's rainfall characteristics.Nevertheless, the localized region of the study does not necessarily reflect the DSD characteristics and optimal QPE relations of the other regions of the country, and it will be very interesting to analyze the intercomparison between the archipelago's other radars representing various regions of the country.

Fig. 1 .
Fig. 1.(a) Radar network in the Philippines; and (b) Area of the study and location of instruments used.
Fig. 5. Reflectivity (dBZ) data (a) before and (b) after attenuation correction from sample dated 08:30 11 August 2018 UTC at 0.5° elevation angle.The red circles highlight some visual difference between the two images.

Fig. 6 .
Fig. 6.Comparison between observed rainfall (R in mm) measured by the rain gauge (x-axis) and values calculated from the Parsivel disdrometer (y-axis) and their accumulated rainfall (mm) for the whole wet season, shown in terms of: (a) scatter plot of rainfall values; and (b) accumulated rainfall.

Fig. 8 .
Fig. 8. Histograms of Dm (blue, in mm) and log 10 N w (red, in mm -1 m -3 ) for (a) the wet season and (b) Yagiand the Distribution of log 10 N w (in mm -1 m -3 ) and D m (in mm) for (c) the wet season and (d) Yagi using 1-min Δt.Also shown are the stratiform-convection separation line and the range values of maritime and continental convection from the study ofBringi et al. (2003Bringi et al. ( , 2009)).

Fig. 9 .
Fig. 9. Scatter plots of N w (10 3 mm -1 m -3 ) and R (mm hr -1 ) for (a) wet season and (b) Yagi as well as for D m (mm) and R (mm hr -1 ) for (c) wet season and (d) Yagi using 1-min Δt.The solid line and the equation represent the fitted power-law curve using least-squares method.

Fig. 11 .
Fig. 11.Scatter plots of DSD-derived log 10 R (mm hr -1 ) and DSD-derived Z (dBZ) for (a) the wet season and (b) Yagi at 1-min Δt.The broken lines are the separation thresholds for stratiform and convective classification in terms of R (yellow) and Z (red).The black dot represents the mean value while the black lines represent the standard deviation in terms of Z (vertical) and R (horizontal).

Fig. 13 .Fig. 14 .
Fig. 12.(a) Comparison between the Z RADAR and Z DSD and (b) their scatter plots Also shown in the figure is the rainfall rate (bars) from rain gauge co-located with the disdrometer.

Table 1 .
Summary of the datasets used in this study.

Table 2 .
Drop size distribution class bins of the optical Parsivel disdrometer.

Table 4 .
Specific parameter values used to retrieve polarimetric variables from Parsivel disdrometer.Technical specifications of the Tagaytay C-band radar.
Length of Observation for the Study 22:00 UTC 10 August 2018 to 16:00 UTC 11 August 2018

Table 6 .
General statistics of rainfall rates (mm hr -1 ) derived from DSD.

Table 7 .
Comparison between gamma parameter values of this study and studies from neighboring countries.

Table 9 .
Validation statistics of QPEs from various R(Z) relations.