Quantitative Precipitation Estimates Using Machine Learning Approaches with Operational Dual-Polarization Radar Data

: Traditional radar-based rainfall estimation is typically done by known functional relationships between the rainfall intensity (R) and radar measurables, such as R–Z h , R–(Z h , Z DR ), etc. One of the biggest advantages of machine learning algorithms is the applicability to a non-linear relationship between a dependent variable and independent variables without any predefined relationships. We explored the potential use of two supervised machine learning methods (regression tree and random forest) in rainfall estimation using dual-polarization radar variables. The regression tree does not require normalization and scaling of data; however, this method is quite unstable since each split depends on the parent split. Since the random forest is an ensemble method of regression trees, it has less variability in prediction compared with regression trees, but consumes more com-puter resources. We considered several different configurations for machine learning algorithms with different sets of dependent and independent variables. The random forest model was appro-priately tuned. In the test of variable importance, the specific differential phase (differential reflectivity) was the most important variable to predict the rainfall rate (residual that is the difference between the true rainfall rate and the one estimated from the R–Z relationship). The models were evaluated by 10-fold cross-validation. The best model was the random forest model using a residual with the non-classified training set. The results indicated that the machine learning algorithms outperformed the traditional R–Z relationship. Then, we applied the best machine learning model to an S-band dual-polarization radar (Mt. Myeonbong) and validated the result with ground rain gauges. The results of the application to radar data showed that the estimates of the residuals had spatial variability. The stratiform and weak rain areas had positive residuals while convective areas had negative residuals, indicating that the spatial error structure driven by the R–Z relationship was well captured by the model. The rainfall rates of all pixels over the study area were adjusted with the estimated residuals. The rainfall rates adjusted by residual showed excellent agreement with the rain gauge, especially at high rainfall rates.


Introduction
Quantitative precipitation estimates (QPE) are a major area of interest within the field of dual-polarization radar. With the advent of polarimetric radar, QPE algorithms using dual-polarization radar variables have been developed in recent decades [1][2][3][4][5]. The dualpolarization radar observes the differential reflectivity (ZDR, dB), specific differential phase (KDP, km -1 ), and cross coefficient (ρ HV ), as well as the reflectivity (Z, mm 6 m -3 or dBZ). Polarimetric variables help to overcome several issues in QPE, such as miscalibration of the radar transmitter or receiver, attenuation in precipitation, and partial beam blockage. The use of these polarimetric variables can provide improved QPE [6]. In addi-tion, since various microphysical information, such as the shape, size, and number concentration of raindrops, is provided using the horizontal and vertical polarization information, QPE using dual-polarization variables provides higher accuracy than using the reflectivity factor in horizontal polarization [7][8][9].
A simple form of radar-based QPE can be performed by an empirical relation between Z and the rainfall rate (R). Marshall and Parmer [10] introduced the R-Z relationship (Z = 200 R 1.6 ), which explains the empirical relationship between Z and R; however, the R-Z relationship is sensitive to the variability of the drop size distribution (DSD, N(D) in m -3 mm -1 ), which causes the uncertainty of the QPE using the R-Z relationship [11,12]. To deal with the uncertainty in the R-Z relationship, rainfall estimation based on dualpolarization radar variables that provide various information on raindrops was proposed [13][14][15].
ZDR is a good measurement of the median volume diameter, and KDP is linearly related to the rainfall rate as it has a lower moment than Z. As a result, rainfall estimation using these variables is more robust with respect to the variability of N(D) [13,16]. It is possible to significantly improve the over/underestimation in rainfall estimation using the R-Z relationship in a strong rainfall rate of 10 mm h -1 or more [17]. These empirical relationships have limitations in explaining the complex nonlinearity between R and radar variables, which leads to errors in rainfall estimation.
Recently, researchers have shown an increased interest in QPE based on machine learning using remote sensing data [18][19][20]. Machine learning (ML) methods, such as decision tree (DT), random forest (RF), and artificial neural networks (ANNs), are techniques for discovering the relationship between independent variables and a dependent variable based on sample data without any preliminary assumptions, including linearity. DT is divided into classification tree (CT) for a discrete dependent variable and regression tree (RT) for a continuous dependent variable [21]. DT can take into account interactions and nonlinearity between variables.
RF is an ensemble method that consists of a number of DTs and shows more accurate prediction performance than a single DT [22]. RF can reduce the variance by lowering the correlation between DTs with randomly selected independent variables. RF can be fitted in parallel, as several DTs are independently generated. Ouallouche et al. [20] performed rainfall estimation based on RF using data retrieved from the satellite, such as the cloud top height, cloud top temperature, cloud phase, and textural features.
As a result, the rain rates estimated by RF greatly agreed with those measured by rain gauges. Kusiak et al. [23] applied five data-mining approaches, including RF and DT, to estimate the rainfall using rain gauge data and Z measured from the Doppler WSR-88D radar of the National Weather Service's Next-Generation Radar (NEXRAD) system. They compared the statistics over the methods, but did not compare them with the empirical relationships. The neural network model showed good performance with the lowest mean absolute error (MAE), and the results were lower in order of support vector machine (0.19), k-nearest neighbor (0.22), CT and RT (0.26), and RF (0.27).
ANNs are ML algorithms that have been inspired by the human neuron-synaptic neural network structure. ANNs are actively applied to atmospheric remote sensing data, as this is effective in extracting characteristics and trends of complex data structures and is suitable for modeling non-linear relationships [24]. Chiang et al. [25] estimated rainfall with a recurrent neural network (RNN) using Z measured from the C-band radar for typhoon periods. The RNN produced better hourly rainfall estimates than those from R-Z relationships in terms of the root mean square error (RMSE).
As a result of the comparison of 48-hours rainfall accumulations, the rainfall estimates obtained from the R-Z relationship were underestimated with a relative bias larger than -45%, while those from the RNN had a relative bias within ±5%. Chen et al. [26] proposed a deep neural network (DNN) approach in rainfall estimation using simulated dual-polarization radar variables based on the N(D) measured from disdrometers. The rainfall rates based on the DNNs model were almost consistent with the rainfall rates computed directly from the N(D).
They compared the hourly rainfall estimates based on the proposed algorithm using Colorado State University-Chicago Illinois (CSU-CHILL) radar data with 11 rain gauges and showed excellent agreement between the estimates and the measurements from the gauges. Although ANNs have been popularly used in a variety of applications due to the advantage of describing nonlinearity between variables, the main shortcoming we encounter with these methods is the black box problem, which makes interpretation of the process difficult. For instance, these methods provide little insight into how the independent variables influence the learning and prediction processes. On the other hand, treebased ML methods (DT and RF) provide ease of interpretation with determination of the variable importance.
Little research has been done on ML-based rainfall estimation using dual-polarization radar. Using the dual-polarization variables allows us to incorporate microphysics information, such as the shape and the number concentration of raindrops into the rainfall estimation. The objective of this study is to improve the accuracy of rainfall estimation based on polarimetric radar parameters using machine learning methods-specifically, tree-based methods (DT and RF).
We used observed drop size distributions, N(D), measured using a two-dimensional video disdrometer (2DVD) to simulate rainfall intensity and radar variables. The ML models were trained with these simulated R and radar variables and cross-validated to check the degree of fitting. The best ML model was independently applied into Mt. Myeonbong (MYN) S-band dual-polarization radar data for rainfall estimation. The estimated rainfall rate was verified using the rain gauge data of automatic weather stations (AWSs) within the radar observation range.

Training Dataset: 2DVD data
In this study, 2DVD data were used to train ML models. The 2DVD is an optical instrument that detects precipitation particles and uses two orthogonal cameras to detect the shadow of particles falling into the observation area. Microphysics information, such as the diameter of particles (D, mm), fall velocity (Vf, m s -1 ), and the axis ratios can be obtained by measuring the shadow of precipitation particles [27]. This information can be contaminated by observing particles that fall into the observation area after being hit by a disdrometer and broken, or by the mismatching of particles in the image processing.
These outliers are removed by comparison with the empirical relationship between D and Vf [28]. In addition, an N(D) that has one or more channels with a zero number concentration is considered as discontinuous N(D) and eliminated [29]. A total of 41 diameter channels of N(D) from 0 mm to 10.25 mm at 0.25 intervals were used. A 1-min rain rate (R, mm h -1 ) was calculated from quality-controlled 1-min N(D) using Equation (1): where dD is the diameter interval at each diameter bin. Table 1 shows the observation locations, observation periods, the number of 1-min N(D), and the maximum 1-min rain rate from the 2DVD data. The total number of training data was 51,302, measured in Oklahoma (OKL, USA), Daegu (DAE, ROK), Boseong (BOS, ROK), and Jincheon (JIN, ROK) to secure the diversity of microphysical processes. We used data observed in spring (April 2019) and autumn (October 2018) as well as summer data (May to September) for seasonal variety. The 2DVD data in Oklahoma were obtained by the National Severe Storms Laboratory (NSSL), National Oceanic and Atmospheric Administration (NOAA), and the data in Jincheon were provided by the Weather Radar Center (WRC), Korea Meteorological Administration (KMA). The other data were collected by the Center for Atmospheric REmote sensing (CARE), Kyungpook National University (KNU).
We discarded the time if the radar reflectivity was greater than 55 dBZ in order to exclude hail particles in the analysis [6]. We also removed the time if the rainfall rate was less than 0.1 mm h -1 because the disdrometer typically underestimates the rainfall rate when small drops (diameter < 0.7 mm) are dominant [30]. The maximum rainfall rate was larger in OKL than the sites in ROK. The median rainfall rate (radar reflectivity) varied from 0.99 mm h -1 (22.92 dBZ) to 1.78 mm h -1 (27.72 dBZ). Table 1 shows the clearly different statistical characteristics of rainfall in different climates (OKL vs. ROK, different regions in ROK) [29]. Unlike the maximum rainfall rate, the maximum reflectivity showed less discrepancy. These characteristics certainly have impact on the ML models and will be discussed later.  [31], and 5-min time average values of Zh, ZDR, and KDP were additionally used to consider the movement of the precipitation system (Zh 5min, ZDR 5min, and KDP 5min). The T-matrix method used to calculate the polarimetric radar variables is one of the most widely used tools for computing light scattering by non-spherical particles based on directly solving Maxwell's equations. This approach can simulate theoretical radar measurements for homogeneous and rotationally symmetric non-spherical particles. Backward scattered fields yield ZH, ZDR, and ρ HV , while forward scattered fields produce KDP [17]. The control conditions and the values used in this study are shown in Table 2. The radar wavelength was 11.01 cm (S-band), and the elevation angle of the radar was set at 0°. The raindrop shape formula suggested by Thurai et al. [32] was used.

Operational MYN S-band Dual-polarization Radar Data
The weather radar used in this study was the MYN S-band dual-polarization radar operated by the Korea Meteorological Administration (KMA). Table 3 shows the detailed specifications of the MYN radar. A volume scan of nine elevation angles (0°, 0.39°, 0.83°, 2°, 2.88°, 4.06°, 5.67°, 7.88°, and 10.94°) was performed every 10 min with a beam width of 0.92°. The measured parameters were ZH, ZDR, KDP, ρ HV , etc. ZH and ZDR were calibrated through post-processing. The averaged ZH calibration bias was calculated by the self-consistency principle of ZH and KDP, and the averaged ZDR calibration bias was conducted by comparing the ZDR and ZH radar measurements with the theoretical relationship between the same parameters simulated by the 2DVD [33,34]. The averaged calibration bias of ZH was -5.0 dBZ, and the averaged calibration bias of ZDR was 0.03 dB. Rain gauge data from a total of 192 rain gauges in KMA AWSs within the MYN radar observation range (150 km) were used to verify the radar QPE ( Figure 1). Each AWS was equipped with two sizes of tipping-bucket rain gauge, 0.1 and 0.5 mm, which measure the 1-min R. In this study, the 10-minute average rain rate was used to match the time resolution with the radar data observed at 10-minute intervals, and missing values were excluded when calculating the 10-minute average rain rate. We analyzed six rainfall cases for QPE and verification. The rainfall events included stratiform rain (Cases 1, 2, 3, and 6) and convective rain (Cases 4 and 5) from 2017 to 2018 (Table 4).

Machine Learning
In this study, RT and RF were used for rainfall estimation. Figure 2 shows a schematic diagram of the RT. We defined N to be the number of RTs, and M to be the number of independent variables. A node was divided into sub-nodes with the lowest variance [21]. A recursive binary partition was conducted for each node until a stop condition was met. The most important independent variable was placed at the top of the tree as a root node. The node divided from the root node is called the intermediate node, and the node that reaches the last is called the terminal node. is an independent variable, and denotes an optimal splitting criterion value for the corresponding independent variables. The ̅ denotes the average value of the data belonging to the jth terminal node.
The RF was based on N RTs with N bootstrap samples (see Figure 3). The bootstrap samples were generated by sampling with replacement. Each RT was grown with the splitting rules using the different independent variables selected randomly. The final prediction was given by the average of the predictions from all RTs. In the RF, the importance of independent variables was measured through the increase of the node purity. The independent variable with the highest increase of node purity played a major role in the prediction.
RF can be optimized by tuning two parameters when generating RTs. One is the number of RTs (ntree = N), and the other is the number of independent variables that are randomly sampled (mtry < M). Liaw and Wiener [35] suggested that ntree was 500, and mtry was √ for classification, as opposed to a third of M for regression. Kühnlein et al. [19] compared out-of-bag (OOB) errors by changing the ntree and mtry to improve the predictability of RF, which was used to select important independent variables and to compute the error of the unbiased estimate [22]. To determine the optimal values with the lowest OOB error in this study, we considered the range of values for each tuning parameter, from 400 to 700 for ntree and from 1 to the number of independent variables for mtry. The optimal ntree and mtry were applied to each model.

R-Z Relationship
For validation of the ML models used in this study, we derived the empirical rainfall estimation relationships based on dual-polarization radar variables. Several R-Z relationships were considered with different thresholds of KDP and ZDR. The R-Z relationship calculated using all training data is shown as Equation (2), and Equation (3) was retrieved with the data below the threshold of KDP and ZDR (KDP < 0.04° km -1 and ZDR < 0.3 dB). Equation (4) was constructed with the data above the thresholds of KDP and ZDR (KDP ≥ 0.04° km -1 or ZDR ≥ 0.3 dB). The equations are assumed to have a power-law (Y = aX b ), and the parameters a and b were estimated using the weighted total least squares method [36].

ML-Based Estimation
The 2DVD data above the thresholds of KDP and ZDR (KDP ≥ 0.04° km -1 or ZDR ≥ 0.3 dB) were used as training data. The dual-polarization radar parameters simulated by the Tmatrix and the 5-min time average values of parameters were used as the independent variables. In this study, we investigated the impacts of three factors on the estimation accuracy. First, three types of dependent variable are considered: R2DVD calculated by N(D) (M1), the residual, ε = R(Z h ) -R 2DVD , between R2DVD and R-which was computed from Equation (4) (M2)-and the normalized residual, ε ̅ = ε/R(Z h ) ̅̅̅̅̅̅̅̅ , (M3).
Second, two different groups of the independent variables were used, and the difference between the two groups included KDP. We denote the two groups by KY for the group with KDP and KN for the group without KDP. Lastly, data binning was implemented to group individual observations into specific bins determined by reflectivity (Table 5), allowing us to train the models locally. The local training with the binned training data and global training with the entire training data are denoted by CY and CN, respectively. In the local training (CY) with RF, ntree and mtry with the lowest OBB error for each model were utilized. KDP and ZDR provide more accurate rainfall estimation than R(Zh) for heavier rainfall, whereas the improvement is not often significant for lighter rainfall due to the noises of KDP and ZDR [13,[37][38][39]. Silvestro et al. [40] proposed a rainfall estimation algorithm that makes use of the best empirical relationships depending on the thresholds of KDP and ZDR, which outperformed R(Zh) for real-time applications. In this study, the rainfall estimation based on this algorithm with different thresholds of KDP and ZDR was performed, which is shown in Figure 4. When KDP was greater than 0.04° km -1 , KY was utilized for rainfall estimation, KN was used if KDP was less than 0.04° km -1 , and ZDR was 0.3 dB or more. The rainfall was estimated by the R-Z relationship (Equation (3)) when both KDP and ZDR were less than the thresholds. The models used in this study are summarized in Table 6.  The trained models were verified by 10-fold cross-validation. Six statistics were used to assess the performance of the ML models: The root mean square error (RMSE), mean absolute error (MAE), bias, correlation coefficient (CORR), coefficient of efficiency (COE) [41], and normalized error (1-NE), where N is the number of observations in test data, Rest represents the estimated rainfall rate, and Robs is the observed rainfall rate.

Application to Operational Radar Data
The ML model with the highest accuracy in the 10-fold cross-validation using 2DVD was applied to the rainfall estimation using the operational dual-polarization radar data. The operational radar data used were conducted using the Hybrid Surface Rainfall method (HSR) of the MYN S-band dual-polarization radar data. The HSR is a technique of generating a rainfall field using the data of the lowest elevation angle that is not affected by ground clutter, beam blockage, and non-meteorological echoes. It is applied to the radar data in polar coordinates, and we selected the rain field using the threshold and calibrated for the radar bias (ZH and ZDR).
When determining the rain field, that is, eliminating non-meteorological echoes or removing artifacts, the threshold values of the following texture were used as follows [42,43]: 0.95 for ρ HV , 0.1 for (ρ HV ), 4.0 dB for (ZH), 4.0 dB for (ZDR), and 15.0° for (Φ DP ). Here, ( ) is the radial texture of variable with a window size of 10. The HSR data at polar coordinates are converted to Cartesian coordinates [42].
The ML-based rainfall estimation was performed for each grid point in the Cartesian coordinates. Similar to the training data, the 5 × 5 km spatial average of Zh, ZDR, and KDP were considered as the independent variables by taking into account the movement of the precipitation system. The estimated rainfall rate was verified by the rain gauges in the AWSs within the radar observation range (150 km). The rainfall rate estimated from the radar grid point closest to the AWS was compared with the average rainfall rate measured by the rain gauges. Figure 5 shows the fitted tree models for M1KYCN and M1KNCN (with and without KDP). See Table 6 for the symbols such as M1, KY, CN, etc. We found that KDP plays an important role in rainfall estimation with M1KYCN because it is used as a splitting criterion at the root node. On the other hand, KDP is not shown in the tree based on M1KNCN, while Zh is the most crucial for the model. The increase of node purity in RF for CN is shown in Table 7. To account for potential overfitting in the random forest, we optimized the tuning parameters, the number of RT (ntree), and the number of independent variables randomly sampled (mtry). The values were expressed according to the inclusion of KDP for each model. For M1KY, the node purity rose the most (713,059) when KDP was used as the criterion for split. Then, in the order of Zh (426,349) KDP 5min (311,490), and Zh 5min (147,971), the increase of node purity was higher. In the case of M1KN, the node purity increased the most when the node was divided by Zh. This indicates that Zh played the most major role in the estimation of the rainfall rate.

Rainfall Estimation from Simulated Dual-Polarization Variables
Unlike M1KY, the increase in the node purity of ZDR was the highest in the M2KY. This indicates that the errors, which cannot be explained by the R-Z relationship, were most closely related to ZDR. The second most important variable in M2KY is ZDR 5min. Similar to KY, ZDR was the most important variable in M2KN, and the importance was the highest in the order of Zh and ZDR 5min. The tendency of increasing node purity of M3 was similar to that of M2 regardless of whether KDP was included or not. As expected, ρ HV was mostly shown to be less important, as it has a very small fluctuation in rainfall cases.  Table 8 presents the increase in the node purity of CY. In M1KY, KDP is shown as the most important variable in all reflectivity intervals as in CN (Table 6), and Zh is indicated as the second most important variable. Similar to CN, M2KY and M2KN always show ZDR as the most important variable. The result of M3 is omitted since it appears with the same tendency as that of M2. The scatterplots of the 10-fold cross-validation for the RT are displayed in Figure 6. The discontinuity of the estimated value is a prominent problem of RT. For the weak R2DVD (R2DVD < 20 mm h -1 ), there is a continuity of data because rainfall is often estimated using the R-Z relationship when the R2DVD is weak and the ZDR and KDP are small. In M1CY (Figure 6b), since the training data are divided by the reflectivity interval, the discontinuity problem is rectified, and the 1-NE value increased 7.08% compared to M1CN. M2 (Figure 6c,d) and M3 (Figure 6e,f) demonstrate similar results and showed a more continuous value than M1, with a higher positive CORR value (CORR > 0.98). Additionally, CY appears lower in the RMSE, MAE, and bias values compared with those of CN, and the CORR, COE, and 1-NE values are higher.   Figure 6. In M1CN (Figure 7a), the CORR value and the COE value are higher than 0.98, even though it pre-sented the worst statistics among the RF models. Compared to M1CN, the underestimation in the strong R2DVD (R2DVD > 80 mm h -1 ) was corrected in M2CN and M3CN ( Figure  7c,e), and the RMSE reduced by 4.93%. The RMSE, MAE, and bias decreased in M1CY (Figure 7b) compared with in M1CN, whereas the RMSE, MAE, and bias increased in M2CY and M3CY (Figure 7d,f).
M2CN outperformed the other models with an RMSE value of 0.598, MAE value of 0.255, CORR value of 0.995, and 1-NE value of 90.99%. The M3CN shows a similar performance with M2CN. There was no significant improvement between CN (left panels) and CY (right panels). Therefore, the ensemble effect by CY was less affected in the RF-based models compared with in the RT-based model. This is because RF itself is based on the ensemble. The interesting thing is that CN in RF even showed a slightly better score than CY. This could be explained by the fact that CN uses more training data and is free from possible local features in CY. The scatterplots of R estimated by the empirical R-Z relationship are shown in Figure  8. The rainfall rate estimated by R(Zh) calculated from the entire training data (Figure 8a) presents a large dispersion overall and a tendency to underestimate at a rainfall rate of 60 mm h -1 or higher. When the rainfall rate was estimated based on Equations (3) and (4) according to the threshold of KDP and ZDR (Figure 8b), the RMSE and MAE increased compared to the case of using one R(Zh) due to overestimating at 40-100 mm h -1 ; however, the underestimation of the high rainfall rate (R > 100 mm h -1 ) was resolved.
Estimation using KDP (Figure 8c) had the lowest RMSE value (1.509) and the highest CORR value (0.974) and COE value (0.939) among the estimations using empirical relationships; however, there was still underestimation overall. On the other hand, in the case of the relationship using Zh and ZDR (Figure 8d), the highest 1-NE value (85.24%) and the lowest MAE value (0.418) are presented, but most of the R2DVD values were overestimated. Among the ML models, the model that showed the lowest performance was M1CN-RT, with an RMSE value of 1.700 and CORR value of 0.961 (see Figure 6). This indicates that all ML-based models outperformed these empirical relationship-based approaches.

Rainfall Estimation from Operational Radar
In the results of the 10-fold cross-validation of the ML models, the M2CN-RF and M3CN-RF showed the best and similar performance in terms of the RMSE and CORR. In this section, we chose M2CN-RF for rainfall estimation with MYN S-band dual-polarization radar data. The HSR images of the rainfall rate of Case 1 and Case 2 are displayed in Figures 9 and 10. R(Zh) was estimated by Equation (3) when KDP and ZDR were below the threshold value and by Equation (4) when they were above the threshold value (Figures  9a,10a).
The M2CN-RF (Figures 9c,10c) adjusted the rainfall rate from R(Zh) by applying the residuals. In the HSR rainfall image, the gray region (ε = 0) in Figures 9c and 10c correspond to the area in which R(Zh) was replaced due to the threshold values of KDP and ZDR. A positive ε value indicates that R(Zh) was overestimated, while a negative ε value indicates that R(Zh) was underestimated. For reference, the HSR of the radar observed variables is shown in Figures 9(d)-(f) and 10(d)-(f).
In Case 1 (1000 LST on August 14, 2017), which is a stratiform case, the area of ε = 0 is wider due to lower values of KDP and ZDR. Although the R(Zh) field is highly correlated with the Zh field, the ε field is less correlated with the Zh field. The large (smaller) value of ZDR is related to the larger positive (negative) ε. The ε is less correlated with KDP in this stratiform case since most of the KDP value is smaller. In Case 2 (0730 LST on September 11, 2017), which is a stratiform rain event with embedded convection, the area of ε = 0 is smaller than that in Case 1. Similar to Case 1, the positive ε area is highly correlated with the higher value of ZDR. Interestingly, the region of higher KDP with Zh > 50 dBZ in the south direction had the largest negative values of ε.
In both cases, it can be seen that ε appeared in space with significant structure, indicating that the error generated from the R-Z relationship had a spatial structure. This result is not surprising because inhomogeneous microphysical processes cause the natural spatial DSD variation and result to the spatial structure of residual from R(Z) [11,12]. The ML model (M2CN-RF) uses the simulated dual-polarization variables, which do not have instrumental noises of the radar, such as sampling error, beam broadening, beam blockage, and miscalibration of the radar [12].
As a result, we can expect the spatial structure of residual will be masked if the order of the magnitude of instrumental noises is comparable to that of natural variation of DSDs. The random variability of instrumental noise can be reduced by averaging the samples. Since the M2CN-RF is an ensemble-based model, the spatial structure of residual can be clearer as the instrumental noises are diminished by averaging the prediction of each RT. Figures 11 and 12 show the verification of the period average rainfall rate (12 h for Case 1, and 11 h for Case 2) for estimation with the empirical relationships and ML. In Case 1, the RMSE value of ML (R(Zh)− ε RF ) is 1.207 and the CORR is 0.773 (Figure 11a), which shows improved performance compared to those estimated by the empirical relationships. When calculated by Equation (2) (Figure 11b), there is a tendency to underestimate when R > 5 mm h -1 . This underestimation is slightly improved by applying Equations (3) and (4) based on the threshold values of KDP and ZDR, and the CORR value increases (Figure 11c).
Underestimation still appears with the heavier rainfall rate. This underestimation of the rainfall rate estimated by R(Zh) is improved by ML in Figure 11a, adjusting the rainfall rate with ε RF . The estimation based on Equation (5) results in substantial underestimation with a low CORR value and the worst performance (Figure 11d). The results estimated by Equation (6) show a severe overestimation of strong RAWS (RAWS > 5 mm h -1 ) (Figure 11e).   In Case 2, a similar trend as seen in Case 1 appears. R(Zh)− ε RF also outperforms the empirical relationships in all statistics. The underestimations of R(Zh) that appear in RAWS stronger than 15 mm h -1 (Figure 12b,c) are corrected by ε RF (Figure 12a), leading to the decreases of RMSE from 2.59 and 1.89 to 1.74, respectively. Analogous to Case 1, Figure  12d presents an overall underestimation, and the RMSE value also shows the largest value. R(Zh, ZDR) overestimated the strong rainfall rate (RAWS > 25 mm h -1 ) (Figure 12e).
A total of six rainfall events from 2017 to 2018 were used to verify the five different models ( Figure 13 and Table 9). The scatter plots of the event averaged rainfall rate are displayed, and the same color indicates the same event ( Figure 6). The statistics are shown for rainfall types and models (Table 1) On the other hand, R(KDP) tended to underestimate (Figure 13d), and R(Zh, ZDR) overestimated in weak rain with RAWS < 17.5 mm h -1 . According to the rainfall types, R(Zh)− had a higher CORR (0.956), COE (0.905), and 1-NE (77.18%) in the stratiform rain, whereas it had a lower RMSE (0.612), MAE (0.330), and bias (1.041) in the convective rain. The 1-NE of the stratiform (convective) rain varied from 60% (10%) to 77% (56%). The poor performance of the convective rain was due to the smaller size of rain cells. Most of the convective rain showed smaller areas of precipitation and short-lived cells.   Table 9. Accuracy of the rainfall estimation with the different models for stratiform, convective, and all cases. The highest values of statistics are highlighted in bold, and the values in parentheses represent the improvement percentage from the best performance of the empirical relationship. Root mean square error (RMSE), mean absolute error (MAE), bias, correlation coefficient (CORR), coefficient of efficiency (COE) [41], and normalized error (1-NE).

Conclusions
The ML-based rainfall estimation using dual-polarization radar variables was explored with simulated and observed variables. The ML methods, RT and RF, used in this study allowed us to model the nonlinear relationship between the dependent and the independent variables and to identify important independent variables. The ML methods were first trained with the DSDs observed from 2DVD. In this study, we also considered three types of dependent variables (R: M1, the residual ε = R(Z h ) -R 2DVD : M2, and the normalized residual ε ̅ = ε/R(Z h ) ̅̅̅̅̅̅̅̅ : M3), two groups of independent variables (with KDP: KY and without KDP: KN), and two types of training data (categorized with intervals of Zh: CY, and overall data without categorization: CN).
In the CY models of RF, the number of RTs and independent variables used was optimized. As a result of ML using DSDs from 2DVD, the KDP was identified as the most important variable for rainfall estimation in both the M1KY-RT and M1KY-RF, while Zh served as the most significant variable in the M1KN. This is an outcome of the fact that the KDP can be approximated with the closest moment of DSDs to R, that is 4.2-5.6th moment and the Zh with the sixth moment of DSDs [12,14]. In M2 and M3, ZDR was the most crucial to explain the error (or residual) occurring in the R-Z relationships.
The ML methods were compared with the empirical relationships through 10-fold cross-validation. In the cross-validation, R(Zh), KN, and KY were first determined with the threshold values of KDP and ZDR due to the noises of KDP and ZDR for light rain ( Figure  4), and the other models were then subsequently applied. Since the estimation in RT took the average of the terminal nodes, discontinuity in the estimation and underestimation (overestimation) at the strong (weak) rainfall rate within the node was often shown; however, residual ε corrects these problems. Similarly, CY is also improved with a lower value of RMSE. Since RF takes an ensemble of RT, the discontinuity problem in RT disappears; however, underestimation was still found above 100 mm h -1 in the M1CN-RF. Adjusting the rainfall rate with the estimated ε resolved the underestimation issue, and the estimation was close to the true value.
Compared to the empirical relationships, all the ML models showed improved evaluation statistics compared with the R-Z relationships. Even the worst ML model (M1CN-RT) showed meaningful improvement of the RMSE value (1.700) compared with the R-Z relationships (RMSE of 2.125 and 3.009). The M2CN-RF outperformed all the empirical relationships and presented a higher CORR value and lower RMSE value. The M2CN-RF, the model with the best performance among the ML models trained with 2DVD data, was applied to the MYN S-band dual-polarization radar data.
In the stratiform case (Case 1), most of the ε values were zero in the weak rainfall rate region, while the ε values were positive in the region of larger ZDR. The negative ε values were large in the convection region, in particular, in the region of larger KDP values in CASE 2. In addition, when estimated by the R-Z relationship, a significant underestimation was shown in heavier rainfall regions and was corrected by ε that was estimated by the ML models. In addition, the significant spatial structure of ε appeared and was highly correlated with ZDR positively, and KDP negatively. The evaluation with six rainfall events indicated that the ML model outperformed the empirical relationships regardless of the rainfall type (i.e., stratiform or convective). The statistics according to the rainfall type show that the ML-based QPE for stratiform cases had a higher CORR, COE, and 1-NE compared with the convective cases.
There was a dependency of the estimation accuracy on the trained data set. When we trained the RF model with DSD data from the DAE that was nearest to the MYN radar, the higher rainfall rate was systematically underestimated, likely due to the limited range of rainfall intensity. A recent study showed a significant discrepancy of DSD characteristics due to different climate and main forcings [29]. The addition of the OKL data significantly extended the rainfall prediction range, particularly in higher ranges, and resolved the underestimation of higher rainfall. As a result, we added additional DSDs data from OKL, BOS, and JIN in the training data set to broaden the variability of DSDs. This improved the accuracy in the overall statistics.
Through this study, we investigated the potential application of rainfall estimation based on ML using polarimetric radar data. We envision that more accurate rainfall estimation can be achieved by applying the RF model considering ε trained with 2DVD data to a large number of operational radars. In particular, the ML model improved the estimates in the heavy rain region, which were underestimated in the empirical relationship. This approach would be useful in the analysis and forecasting of severe weather. Future studies could include extending ML models to various radars and rainfall cases in different weather conditions.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.