Geostatistical Based Models for the Spatial Adjustment of Radar Rainfall Data in Typhoon Events at a High-Elevation River Watershed

: Geographical constraints limit the number and placement of gauges, especially in mountainous regions, so that rainfall values over the ungauged regions are generally estimated through spatial interpolation. However, spatial interpolation easily misses the representation of the overall rainfall distribution due to undersampling if the number of stations is insufficient. In this study, two algorithms based on the multivariate regression-kriging (RK) and merging spatial interpolation techniques were developed to adjust rain fields from unreliable radar estimates using gauge observations as target values for the high-elevation Chenyulan River watershed in Taiwan. The developed geostatistical models were applied to the events of five moderate to high magnitude typhoons, namely Kalmaegi, Morakot, Fungwong, Sinlaku, and Fanapi, that struck Taiwan in the past 12 years, such that the QPESUMS’ (quantitative precipitation estimation and segregation using multiple sensors) radar rainfall data could be reasonably corrected with accuracy, especially when the sampling conditions were inadequate. The interpolated rainfall values by the RK and merging techniques were cross validated with the gauge measurements and compared to the interpolated results from the ordinary kriging (OK) method. The comparisons and performance evaluations were carried out and analyzed from three different aspects (error analysis, hyetographs, and data scattering plots along the 45-degree reference line). Based on the results, it was clearly shown that both of the RK and merging methods could effectively produce reliable rainfall data covering the study watershed. Both approaches could improve the event rainfall values, with the root-mean-square error (RMSE) reduced by up to roughly 30% to 40% at locations inside the watershed. The averaged coefficient of efficiency (CE) from the adjusted rainfall data could also be improved to the level of 0.84 or above. It was concluded that the original QPESUMS rainfall data through the process of RK or merging spatial interpolations could be corrected with better accuracy for most stations tested. According to the error analysis, relatively, the RK procedure, when applied to the five typhoon events, consistently made better adjustments on the original radar rainfall data than the merging method did for fitting to the gauge data. In addition, the RK and merging methods were demonstrated to outperform the univariate OK method for correcting the radar data, especially for the locations with the issues of having inadequate numbers of gauge stations around them or distant from each other.


Introduction
Hydrologic predictions in terms of surface runoff are critically important to disaster mitigation as well as water resource planning and management. Tremendous efforts have been dedicated to the development of numerical models by researchers over the years in attempts to accurately simulate the hydrologic responses of targeted areas during severe weather events. With the increase in computational power and availability of high resolution spatial and temporal data, physically based distributed models have rapidly gained popularity over empirically based lumped models. Complex governing equations can now be solved using advanced numerical methods without compromising the computational efficiency. Meanwhile, an increase in streamflow prediction accuracy can be achieved with the inclusion of the detailed geophysical and meteorological characteristics of areas of interest. The performances of distributed and lumped models have been assessed (e.g., [1][2][3][4]). Studies on the implementation of the distributed model for a variety of hydrologic applications can be found in [5][6][7][8][9][10].
Among all the required inputs for runoff modeling, rainfall data undoubtedly have the most critical influence on the prediction results. In order to fulfill the requirement of high-resolution input data for distributed simulations, dense gauge networks must be built to observe the actual spatial distribution of rainfall. However, establishing dense gauge networks, especially in mountainous regions, is impractical due to the high cost of instrumentation and required routine maintenance services [11,12]. In addition, geographical constraints such as highly sloped terrain also limit the number and placement of gauges within a network. Rainfall values over the ungauged regions are generally estimated through spatial interpolation. As indicated by Duncan et al. [13], the interpolated results were heavily dependent on rain gauge network density. Regardless which interpolation technique is chosen, the technique can easily misrepresent the overall rainfall distribution due to undersampling if the number of stations is insufficient. Sampling errors produced from rain gauge networks with various densities have also been investigated by other researchers [13][14][15][16].
With the recent advancement in observation technology, remote sensing has become an invaluable alternative to field sample collection because of its ability to economically acquire the detailed spatial distribution of target variables over large areas in a timely fashion [17][18][19]. Differing from satellite-based precipitation measurements [20][21][22], weather radar is considered to be one of such systems widely used to provide multiple types of meteorological products (e.g., precipitation, hail index, wind profile) [23]. These products are typically derived from reflectivity measurements and stored in gridded formats and short intervals in time, which are particularly useful for real-time disaster warning operations. Radar precipitation has had a long history of being used in flood forecasting and other types of hydrologic predictions [1,11,[24][25][26][27][28][29][30][31]. While the precipitation estimates provide the spatial and temporal details of a rainfall event, it is well known that the data are subject to uncertainties and systematic errors [11,25,[32][33][34][35][36][37]. Hence, proper adjustments are required prior to their use as inputs for any simulation tasks. As pointed out by Wilson and Brandes [36], the main source of error for radar rainfall occurs during the conversion process. The precipitation estimates are converted from the reflectivity measurements using the Z−R relationship, where Z is the radar reflectivity and R is the radar rainfall rate, under the assumption that the rain drop distribution is exponential and the vertical air flow is relatively small compared to the terminal velocities of rain drops. In reality, the drop size distribution and air motion will most likely not be as assumed and vary continuously throughout an event. It would not be feasible to adjust the coefficients of the Z−R relationship for every event.
As rainfall recorded by rain gauges is commonly considered as ground truth, a great number of methodologies that combine both gauge observations and radar estimates to generate rain fields have been proposed. In principle, by utilizing data from both sources, the accuracy of point measurement and high-resolution spatial information can be integrated to potentially obtain a measurement more accurate to the actual rainfall distribution over targeted areas for the use of hydrological modeling and flood forecasting [5,30,38]. Improvements to the accuracy of rainfall estimations and hydrological modeling with rainfall inputs through the combination of radar and rain gauge observations have been reported by Kim et al. [39], Gourley and Vieux [40], and others. Studies, including Morin et al. [41], Rosenfeld et al. [42], and Smith et al. [43], also addressed the rainfall accuracy related to the adjustment of radar data using gauge measurements and affecting parameters. Researchers and governmental agencies have adopted that the mean field bias correction method is a quick and simple way to calibrate radar rainfall in real time operations [24,26,30,36,44,45]. The ratios of gauge to radar values at the locations of all gauges in a network are averaged, and the averaged correction factor is then multiplied uniformly to the entire radar field. Considering the mean field bias between radar areal estimates and the point measurements from gauges, Kalman Filter's procedure has been applied to adjust the NEXRAD (next generation weather radar) radar rainfall forecasting system [26,46]. Focusing on removing the systematic errors or bias, Rosenfeld et al. [42] proposed a probability matching method for radar adjustment. Cole and More [1] used an integrated multiquadric estimation technique [47] for the dynamic generation of gridded gauge-adjusted rainfall fields.
Geostatistically based approaches in recent years have become increasingly preferred for the interpolation of rainfall. Kriging or ordinary kriging (OK) is one of such approaches that has been considered the best linear unbiased estimator [48]. Some commonly seen forms of kriging that have shown excellent rainfall prediction results using both gauge precipitation and radar estimates or other secondary variables include simple kriging with varying local means (SKlm) [49], kriging with external drift (KED), regression-kriging (RK) [50,51], universal kriging (UK), and co-kriging (CK). SKlm, KED, RK, and UK are very similar in terms of their mathematical formulations. The primary differences between these multivariate techniques are the steps that are taken to solve the kriging weights. Compared to the others, the ordinary CK is a more complex and computationally demanding extension of kriging. Not only the does the correlation of each of the primary and secondary variables have to be analyzed, the cross correlation between the residuals of the variables must also be examined [12]. The various forms of multivariate kriging for rainfall interpolation in distributed hydrologic modeling have been assessed in the past. Goovaerts [52] explored the use of the digital elevation model (DEM) as an external predictor for spatial rainfall prediction. The results showed that SKlm, KED, and CK outperformed the conventional Thiessen polygon [53], inverse distance weighting (IDW), OK, and linear regression. Haberlandt [54] included both radar estimates and elevation values as external predictors for rain field interpolation. The results also confirmed that multivariate methods clearly outperformed univariate techniques. The comparison study for the KED and RK methods were presented by Hengl et al. [55]. Zhang and Srinivasah [45] employed the RK and SKlm approaches for validating NEXRAD data. Later, Schuurmans et al. [56] studied the effect of the extent of the interpolation domain on the spatial rainfall prediction accuracy. The performances of KED and CK using radar precipitation as a secondary predictor were compared against OK. Ehrat [38] introduced a geostatistical merging method for spatial rainfall estimation using both radar and rain-gauge data. Other tested case studies can also be found in [57][58][59][60]. However, the above mentioned geostatistically based methods have not yet been tested and examined for their performances applying to a high-elevation study area.
Every summer and fall, the frequently occurring regional thunderstorms, seasonal rain (named Meiyu), or typhoons are known to routinely bring large amounts of rainwater in a short period of time to Taiwan. Typhoons invade Taiwan every year and have been observed to cause severe floodrelated property damage and the loss of lives regularly. The study area, Chenyulan River watershed, is located in central Taiwan, in which watersheds are highly prone to flash flooding during the typhoon season, especially after the Chi-Chi earthquake in 1999 [61]. The early warning system that is currently in operation for flood forecasting relies solely on its rain gauge network as the rainfall data source for hydrologic and hydraulic simulations. Due to the rough terrains, high elevation, and heavy forest cover, establishing and maintaining a functional high-density rain gauge network in the study area has been proven to be a challenge. The ultimate goals of this study would be to integrate the developed algorithm into the system and enhance the accuracy of rainfall input with the muchdetailed high resolution space-time rainfall data. The algorithms developed in this study were based on two spatial interpolation-based methods, regression-kriging (RK) and merging, and written in the R statistical programming language [62]. The performance of the method was evaluated on five historical typhoon events (Kalmaegi, Fungwong, Sinlaku, Morakot, and Fanapi). Error analyses of the interpolated results were performed. The corrected radar rainfall values along with the results obtained using ordinary kriging (OK) were compared in their accuracies against the observed rainfall values at the gauge locations. The time series plots of the adjusted rainfall data showing the improved variation trend following the gauge measurements are also presented to confirm the performance of the developed spatial interpolation models.

Study Area and Rainfall Data
The Chenyulan River watershed is located in central Taiwan and is elongated in the north-south direction and encompasses approximately 450 km 2 of land. The area rises to over 3700 meters near Mount Yu and descends northward to an elevation of 400 meters at its outlet. With the central ridge of the island of Taiwan to the east, the averaged elevation in the eastern and southeastern regions of the watershed is above 2500 meters. The location and elevation distribution of the watershed are shown in Figure 1.
The watershed consists mostly of mountainous terrain. The slope of the majority of the land surface ranges from 20 degrees to 60 degrees. The main river, the Chenyulan River, originates from the north ridge of Mount Yu and flows an estimated distance of 42.4 kilometers through a deeply cut valley. Due to a significant elevation drop of more than 2000 meters between the north and south ends of the watershed, the averaged channel bottom slope of the Chenyulan River is nearly 6.75 percent. The overall riverine system is depicted in Figure 1. The combination of highly sloped river channels and steep terrain accelerates the rainfall-runoff transformation, which makes the watershed highly susceptible to flash flooding [63]. The nearly 15% flatter areas are generally near the bottom of the river valley and are where most of the human activities take place (e.g., residential, agricultural, and other types). However, more often than not these areas are also located in the floodplains, such that time for emergency response during severe weather events is greatly limited.
This study incorporates the true gauge observations and raw radar estimates to perform the adjustment of radar rainfall data using geostatistical-based interpolation procedures. The importance of adequate rainfall sampling for rain field interpolation during storm events can never be overemphasized. As aforementioned, accurate rainfall data can enhance the modeling results. In order to achieve adequate sample collection, sufficient and evenly distributed sampling locations are generally required. However, various geographical constraints often dictate the number and location of the sampling sites when designing a rain gauge network. The collected radar rainfall data can serve as a practical solution, providing a wider and more detailed distribution of the rainfall field, but using an effective interpolation tool to correct the rainfall values is required.
The rain gauge network used for this study consists of 27 strategically selected stations, among which 23 are managed by the Central Weather Bureau (CWB) of Taiwan. The remaining four stations are operated by the Water Resources Agency (WRA). The gauge network is shown in Figure 2. The vertical elevation for the gauge stations ranges from 235 to 3845 m, with an average elevation of 1504 m (see elevation data in Table 1 and Table 2). Among the 27 stations, 12 stations are located inside the Chenyulan River watershed and the other 15 gauge locations are stationed in the regions either outside or on the watershed boundary. Additionally, 17 out of the 27 stations have the elevations above 1000 m. Due to the rough terrain restriction, most of the rain gauges in the study area were placed along river channels where easier access to the sites could be gained. In high elevation regions such as the east part of the watershed, establishing gauging stations was almost impossible. Therefore, undersampling in certain areas could be expected.
In regard to rainfall recorded at each gauge station, the CWB preprocesses and cumulates the measured rain amount over regular time intervals (10 minutes and 1 hour) for their managed 23 stations. However, the rainfall data collected in the 4 WRA stations were in an irregular time series format. Extra interpolation efforts were carried out to convert the rainfall measurements in WRA stations into a regular time series format that was consistent with the data from the CWB, including time interval and time stamp format. Then, the regular time series of the rainfall data from all gauge stations could be used directly by the rain field interpolation algorithm developed in this study.  The CWB, WRA, Soil and Water Conservation Bureau (SWCB) of Taiwan, and US National Oceanic and Atmospheric Administration's (NOAA's) National Severe Storms Laboratory (NSSL) had collaborated in the development of the quantitative precipitation estimation and segregation using multiple sensors (QPESUMS) system. Since its deployment in 2002, the system has been used to produce a variety of weather products by integrating radar data with other meteorological observations such as wind speed, lightning strikes, and rain gauge measurement. Through the CWB's QPESUMS system, the radar rainfall data were stored in a gridded format with a bin dimension of 0.0125 degree × 0.0125 degree (latitude/longitude), i.e., a cell size of approximately 1.25 × 1.25 km, and the recording interval was ten minutes. Figure 3 presents an example plot showing a cropped time slice of the unadjusted radar rainfall of Typhoon Morakot. The detailed spatial distribution of the rainfall at a selected instant can be displayed.

Regression Kriging
The raw radar data extracted from QPESUMS contain the rainfall values of the covered area. They were converted from the measured radar reflectivity by applying the Z-R relationship. According to the CWB of Taiwan, the Z-R equation accepted for rainfall conversion in Taiwan is given as where Z is the radar reflectivity in dBZ and R is the radar rainfall rate in mm/hr. One of the selected methods for correcting the R values using gauge rainfall measurements is the geostatistically based regression-kriging (RK) technique, where the spatial variation in a target variable can be presented as the summation of its deterministic (trend) and stochastic (disturbance or residual) components. Regression-kriging is a combination of the generalized linear models with kriging. In RK, the deterministic part of estimation of a target variable is modelled by the linear regression approach using the formulated relationship between the target and auxiliary variables. The stochastic part of the calculation is estimated by kriging based on the regression residuals. In this study, the target variable was the rainfall rate from the gauge stations, and the auxiliary variables were the radar rainfall estimates from the QPESUMS system and the elevation values pertaining to the watershed and its surrounding areas. According to Hengl [64], the variation trends of the measured rainfall rate ( ) at the gauge stations can be commonly related to the affecting parameters using linear regression, as Using the pairs of ( ) and ( ) values and the elevation data, the coefficients , , and are determined by the procedure of generalized least squares regression. With the known regression coefficients, the trend of rainfall variation for a given instant at any ungauged location = , can be obtained using the corresponding radar and elevation data at . A sample plot showing an interpolated trend surface of the rainfall variation is presented in Figure 4a.
In addition to a trend in the rainfall distribution surface, a residual surface has to be produced.
The residual at any ungauged location ( ) is estimated by the kriging technique. The estimated residual surface is a linear combination of the residuals at neighboring gauge stations multiplied by the appropriate kriging weights as described by ̂ = ∑ 1 ( ), where are the kriging weights determined at with respect to . Figure 4b illustrates an example plot of the interpolated residual rainfall surface. A semivariance analysis with fitted semivariogram functions can be performed to determine the weighting coefficients. The semivariance (SV) of the residual is computed using the definition shown below: Here, ( ) = the magnitude of the regional variable (or the residual); ( + ℎ) = the magnitude of the regional variable that is away from ( ) by a distance of h; (ℎ) = SV function; and (ℎ) = the number of pairs of residual variables separated by h. A semivariogram can be generated by plotting semivariances versus distances between ordered data. Figure 5 shows an example of a semivariogram. It should be noted that strong similarities exist between the residuals at shorter distances, as the semivariance decreases with the decreasing distance between gauge stations.  With an established semivariogram, the kriging weights are determined by solving the following system of equations.
where = values obtained from the residual related semivariogram according to the distance between two points denoted by i and j; = values of semivariance based on the distance between an ungauged point p and gauge point i; and μ = Lagrange multiplier. The final adjusted radar rainfall rate, , at any ungauged location can then be determined by combining the calibrated variation trend and, with the use of the kriging technique, the geostatistically estimated residual. We have = + + + ∑ ( ) .

Merging Method
Apart from the RK method, the other approach considered in this study was the merging method, which also jointly uses radar and rain-gauge data. The idea of the merging method is to merge the interpolated gauge data at an ungauged station with the deviations in the observed and interpolated radar values at that station. The effect of reliable spatial distribution from the observed radar rainfall field can be included in the correction process. Again, assuming to be the adjusted radar rainfall rate at any ungauged location , the merging method with the application of kriging technique gives where w are the weighting coefficients, ( ) again is the rainfall rate measured at any gauge station located at , and ∆ represents the rainfall correction obtained from the information relating to observed and interpolated radar values. Similar to Equation (5), the coefficients w are determined by solving a system of equations with known coefficients obtained from a semivariogram constructed with the semivariances of rain-gauge data between any pair of gauge stations. If the effect of ∆ is neglected, Equation (7) is recovered with the ordinary kriging (OK) method. For the determination of the ∆ values, the radar values at gauge locations and the kriging procedure are first used to estimate the spatially interpolated radar rainfall rates at an ungauged location, , using the formula ∑ ( ). Then, the rainfall correction at any ungauged location can be computed as where represents the original radar rainfall rate at the targeted location . ( ), = 1, 2, … , are the original radar rainfall data identified at the gauge locations, , , … , . Again, the kriging weighting coefficients are determined by solving a system of equations similar to Equation (5). However, with the coefficients obtained from a semivariogram, that is assembled with any pair of radar data ( ) at the gauge locations and distance h (See Equation (4)). Because the RK (or merging) method uses the QPESUMS estimates as an auxiliary variable to interpolate the rainfall, the interpolation can be perceived as a process to adjust the QPESUMS data to match the observed rainfall data.

Results
The algorithms, according to the above-described RK and merging methods, were coded in the R language and applied to the Chenyulan river watershed to correct the raw radar rainfall data recorded under five selected severe typhoon events in last decade, which are typhoons Kalmaegi, Morakot, Fungwong, Sinlaku, and Fanapi. The performance of the RK and merging method were evaluated. The corrected results after running the RK and merging codes for the tested typhoon events were cross validated with the gauge observations. In this study, the leave-one-out cross validation, which is widely recognized as one of the most efficient ways to validate a prediction model, was performed to assess the accuracy of the radar rain fields interpolated from the RK and merging models. Their percentages of improvement were compared with those from the usual OK method, which is an interpolation procedure just based on the gauge data. The observed rainfall data from one of the gauge stations were first excluded from the interpolation process to be used as the validation dataset, while the observed data from the remaining stations were utilized with the radar rainfall estimates to interpolate the rainfall at the location of that gauge station. The procedure was repeated until each individual gauge station had been rotated through as a validation location. Finally, the corrected rainfall values were compared to the gauge data for accuracy verification. Three gauge locations-station C1M440 (upstream), station C1I060 (central region), and station C1I080 (downstream)-along the Chenyulan River were selected to detail the comparisons between the time series of the adjusted radar rainfall rates and the observed data.
An error analysis was carried out to show the overall performance of the proposed interpolation techniques. In this study, two indices, the root-mean-square error (RMSE) and the coefficient of efficiency (CE), were computed to quantitatively evaluate the errors produced by the proposed models. The indices were calculated at all gauge locations considering the full event data from the typhoons tested in this study. The CE was calculated as where again represents the position vector of the ith gauge station, ( ) is the jth interpolated (or adjusted) rainfall value at , ( ) is the jth gauge data at , ( ) is the time averaged value of the gauge rainfall data at , and N is the total number of data points. It should be noted that if the interpolation techniques presented in this study made the rainfall adjustments accurately, the resulting CE would be close to unity, as the second term on the right-hand side of Equation (9) would be relatively small. The results of the time series comparisons in the rainfall values, error analysis, and the data scattering plots along the 45-degree reference line, showing the fitness of the corrected rainfall values versus the measured ones, are presented in the following subsections.

Typhoon Kalmaegi
Typhoon Kalmaegi swept through Taiwan in July 2008 (07/16/2008 -07/18/2008). A total depth of nearly 700 mm of rain was recorded within a two-day period. The RMSE and CE values-based on the comparisons between the true gauge data and, separately, the original QPESUMS estimates and those corrected using the RK and merging methods-are summarized in Tables 1 and 2. The RMSE reduction percentage and the CE improvement percentage in reference to those of the QPESUMS data are also listed in Tables 1 and 2 for evaluations of the model's performance. Similarly, the calculated RMSE and CE values, the RMSE reduction percentage, and the CE improvement percentage from the standard OK method are given in Tables 3 and 4 for comparison. A positive reduction percentage in RMSE or a positive improvement in CE indicates that the proposed correction models could produce more accurate rainfall values than the QPESUMS, whereas a negative percentage would mean that the approaches tended to not improve the original QPESUMS data. The results of the error analysis are divided into two groups. The first group (Table 1 and Table  3) includes the gauge stations that are located inside the Chenyulan River watershed where the stations are closer and denser. The stations that are outside (or on) the watershed boundary belong to the second group (Table 2 and Table 4), and are generally distant from each other. The impact of the gauge locations and sampling quality on the error reduction for the adjusted radar rainfall values is shown, as the averaged RMSE reduction percentage for the first group is much higher than that of the second group. For the first group (Table 1), the RK method yields an averaged RMSE reduction percentage of 26.93%, followed closely by a 26.39% reduction with the merging method. Referring to the maximum reduction percentage, both methods can improve the rainfall values, with the RMSE reduced by up to nearly 50%. From Table 3, it can be noted that the averaged reduction percentage of the RMSE from the standard OK method is lower, at 22.24%. For the gauge stations outside or on the watershed boundary, as shown in Table 2, the averaged RMSE reductions are 11.10% and 9.32% for the RK and merging methods, respectively. It is believed that inadequate sampling values were recorded at station C1I100, so the averaged percentage values for the RMSE were calculated by excluding the values from that station. When the OK procedure was used, the averaged RMSE reduction percentage was -9.58% (or -25.03% if the negative percentage value at station C1I100 was included), indicating the interpolated rainfall values are less accurate than the QPESUMS data. This may be a result of insufficient gauge stations in the region outside the watershed boundary limiting the performance of the OK method. The RK and merging models are shown to give more improved radar rainfall estimates than the OK approach does, especially for locations in a region with a less dense gauge distribution. The results of the coefficient of efficiency (CE) are summarized in Tables 1 to 4 present further evidence of the improved performance of the RK and merging methods. As mentioned above, a CE value approaching unity would indicate that the interpolated values are close to the observed gauge values. From Table 1, it is noted that the averaged CE of the QPESUMS data from the gauge stations located within the watershed boundary is about 0.72. With the procedures of data adjustment carried out by the RK and merging models, the averaged CE of the corrected rainfall data by either method improves equivalently to 0.84. The maximum CE value can reach up to 0.95. In comparison, the averaged CE using OK is 0.82 (See Table 3). It can be seen that more accurate rainfall adjustments can be made with the proposed geostatistical methods. In the event of Typhoon Kalmaegi, the averaged CE improvement percentage for stations within the watershed using the approaches of RK and merging are respectively 16.72% and 16.35%. The analysis has shown again that RK outperforms the merging method slightly and that both the RK and merging models are notably more effective than the OK procedure in correcting the unreliable radar rainfall data. Similar to the RMSE-based analysis, as seen in Table 2, the RK and merging methods again generate equally promising results, with averaged CE values of 0.78 for the gauge stations outside the watershed. The averaged improvement percentages are respectively 6.99% and 6.15% from the RK and merging results, representing a lower percentage when compared to the results in Table 1. Again, due to the inadequate sampling values found at station C1I100, its CE and percentage values were excluded in the calculation of all CE-related averaged values. From Table 4, the averaged CE from OK for those outside the watershed locations is, however, 0.65 (if the CE value at station C1I100 (i.e. −7.77) is included, the averaged CE is reduced to 0.09), which is lower than the QPESUMS' CE value of 0.73. Also, for those stations (not counting station C1I100) outside the watershed, two stations are shown to receive small negative percentage values using the RK or merging approaches, while 7 out of 14 stations were found to have negative values when the OK method was used. The averaged percentage improvement in CE from OK is -12.94%. Therefore, again, from the comparisons made to the CE values, the performance of RK and merging are found to be substantially better than that of OK.
The time variations in rainfall depth (hyetographs) at station C1I060 showing the adjusted event rainfall values from RK and merging, the original QPESUMS data, and the gauge measurements for Typhoon Kalmaegi are presented in Figure 6. It can be seen in Figure 6 that the QPESUMS data do not follow the measurements closely and tend to underestimate the local rainfall peaks and overestimate the rainfall values in the recession limb of the event. In contrast, the interpolated rainfall values of the proposed RK and merging methods follow the gauge measurements more closely (capturing the peaks and troughs). The gauge data recorded during the entire typhoon event and the corresponding adjusted radar values were plotted as pairs of data points with the 45-degree reference line overlaid for comparison. The plots are shown in Figure 7a,b,c, respectively, for stations C1M440, C1I060, and C1I080. The unadjusted QPESUMS data were also included. It can be concluded that, at each station, the adjusted radar data agree reasonably well with the gauge data, as indicated by the narrower bandwidth of the data points close to the 45-degree reference line. The unadjusted QPESUMS data, however, are shown to scatter to a wider extent (wider bandwidth). The RK and merging methods are again demonstrated to be able to reliably perform the correction procedures on the raw radar data. Based on the results shown in Figure 7a-c, a better agreement between the gauge data and adjusted radar values can be noticed for stations C1I060 and C1M440. For station C1I080, both the RK and merging methods tend to slightly under adjust the raw radar data for rainfall depths greater than 10 mm (interpolated data points from RK and merging fall within the lower region).

Typhoon Morakot
Typhoon Morakot, with a duration from 08/06/2009-08/10/2009, made landfall on August 7, 2009. As a record-breaking event, Typhoon Morakot dumped tremendous amounts of rain on Taiwan-e.g., nearly 3150 mm of rain was recorded at station 467530. The proposed RK and merging algorithms were tested with this historical rainfall event for verification. The results of the error analysis based on the calculated RMSE and CE indices are summarized in Table 5 and Table 6. For the gauge stations inside the watershed, the averaged RMSE reduction percentages of 39.75% and 37.91% (See Table 5) are achieved by the RK and merging methods, respectively, which are better than that of 32.67% (See Table 3) from the OK procedure. When considering the stations outside the watershed (See Tables 6), the averaged RMSE reduction percentages from RK and merging are respectively 17.93% and 11.48%. However, the OK again gives a negative averaged RMSE reduction percentage (i.e. −16.92%) (see Table 4) and 11 out of 15 stations are not shown to have improved rainfall values, indicating that the OK method is not able to improve the accuracy of the raw radar rainfall data.
In terms of the averaged CE improvement percentage based on the RK and merging models, they are 36.21% (average CE=0.87) and 34.94% (average CE=0.86) (see Table 5), respectively, for the gauge stations inside the Chenyulan River watershed, and 24.50% and 16.51% (see Table 6) for the gauge stations outside the watershed. However, the averaged CE of the OK is reduced to 0.34 (see Table 4), which is less than the QPESUMS' CE value of 0.56. It is also shown that the OK approach produces a negative value in the averaged percentage improvement in CE (i.e. −45.89%). According to the results of the analysis, RK is shown to be the most accurate in interpolating and improving the rainfall data, followed by merging. Similar to the case of Typhoon Kalmaegi, the performance of the OK method for locations outside of the watershed is not as consistent as the RK and merging methods due to inadequate sampling locations and limitations of the OK methodology. The hyetographs covering the event period of Typhoon Morakot at station C1I060 are shown in Figure 8. The plots include the interpolated results from the RK and merging methods. Observations from the gauges and the unadjusted QPESUMS radar estimates are also plotted for comparison. It can be clearly seen that the unadjusted QPESUMS data show poor representation of the real-time varying rainfall values. The adjustment procedure from RK is shown to be able to substantially improve the radar data to fit closely to the gauge measurements, especially within the period from midnight to 1:00 a.m. on August 9, 2009. To better the comparisons of the performances of the RK and merging methods, all the interpolated rainfall values and unadjusted QPESUMS data for the event of Typhoon Morakot are plotted against the gauge measurements in Figure 9a-c for stations C1M440, C1I060, and C1I080, respectively. The 45-degree reference line represents the perfect agreement. As shown in Figure 9ac, the data points of the rainfall values corrected by the spatial interpolation algorithms of RK and merging scatter around the reference line more closely than those of the QPESUMS data. The plots also support the statement made above that the adjusted radar rainfall values at station C1I060 among the three test stations are shown to have the best agreement with the gauge data.

Typhoons Fungwong, Sinlaku, and Fanapi
The results of the cross validation of the adjusted rainfall values by the RK and merging models and their performances, measured by the percentage reduction in RMSE and percentage improvement in CE for the other three typhoon events-namely, typhoons Fungwong (duration: 07/26/2008-07/29/2008), Sinlaku (duration: 09/11/2008-09/16/2008), and Fanapi (duration: 09/17/2010-09/20/2010)-in general, are very similar to those presented in the cases of typhoons Kalmaegi and Morakot. In these three additional test cases, less extended error analysis results are presented. The summaries of the averaged RMSE reduction percentage and the averaged CE improvement percentage for the RK, merging, and OK methods for the gauge stations located within the Chenyulan River watershed for these three typhoon events are provided in Table 7. It should be pointed out that unusually large sampling errors were identified at gauge station C1M440 for Typhoon Fungwong and station C1I070 for Typhoon Fanapi. Thus, the error analysis data from those two stations were not used in the calculation of the averaged values shown in Table 7.  Based on the results indicated in Table 7, the performances of the three interpolation techniques can be ranked from best to worst, respectively, as RK, merging, and OK. It is noticed that the merging and OK interpolation methods were not as effective in correcting the rainfall values for Typhoon Fanapi as with Typhoons Fungwong and Sinlaku, as the averaged percentages for the RMSE reduction and CE improvement are shown to be in a lower or negative range for Typhoon Fanapi when compared to the values from other two events. However, even for the event of Typhoon Fanapi, the RK method can still show a reasonable performance in adjusting the radar data. It is interesting to note that, among the five typhoon events analyzed, the merging method has failed to yield a positive averaged CE value for the case of Typhoon Fanapi.
As examples, the comparisons of the rainfall time series for the events of typhoons Fungwong and Sinlaku at station C1I060 are presented in Figure 10 and Figure 11, respectively. The results in general have further confirmed that the RK and merging methods can perform reasonably well in correcting the radar rainfall values for severe typhoon events.

Discussion
As seen in the previous section, the results and their interpreted findings from the comparison plots of the adjusted radar rainfall values against the gauge measurements and the error analysisobtained RMSE and CE values from the developed RK and merging models-applied to the five selected severe events of typhoons Kalmaegi, Morakot, Fungwong, Sinlaku, and Fanapi-were presented to evaluate the levels of improvement made by the RK and merging methods against the unadjusted QPESUMS rainfall data. In this section, the key features in connection with the performance of the RK and merging methods when compared to the conventional OK method and their effectiveness when applied to the high-elevation stations are discussed.
Based on the performance comparisons measured by the RMSE and CE between the adjusted or unadjusted radar rainfall data with the gauge measurements for the typhoon events tested, it is found that the obtained values for the averaged percentage reduction in RMSE and percentage improvement in CE from the RK and merging methods are significantly larger than those from the OK method, especially at stations that are located outside of the Chenyulan River watershed. Therefore, the RK and merging models are considered to outperform the OK approach and are demonstrated to be more effective than the OK procedure in correcting the radar rainfall data. When examining the difference in performance between the RK and merging methods, the results from the error analysis reveal that RK can produce slightly better adjustment in radar rainfall values than the merging method can.
In this study, the developed RK and merging models were applied uniquely to a high-elevation watershed. As indicated in Table 1 Table 2). It would be interesting to examine and discuss the effectiveness of the RK and merging methods performed for the adjustment of radar rainfall data at a high-elevation study area. Here, the stations with elevations higher than 1000 m are defined as high-elevation stations. To make comparisons, the RMSE and CE values in Tables 1 and 2 and Tables 5 and 6 are used to calculate the averaged percentage reduction in RMSE and percentage improvement in CE for the adjusted radar rainfall estimates from RK and merging separately for stations that are above 1000 m and those below 1000 m.
Considering the case of Typhoon Kalmaegi, the overall averaged values for percentage reduction in RMSE and percentage improvement in CE are 29.3% and 16.7% and 27.99% and 15.67%, respectively, when applying the RK and merging methods to the stations that are inside the watershed boundary and above 1000 m in elevation. For the stations that are inside the watershed boundary but below 1000 m in elevation, the overall averaged percentage reduction in RMSE and percentage improvement in CE from the RK and merging results are 24.56% and 16.74% and 24.8% and 17.02%, respectively. It can be seen that, when measured by the statistical values of RMSE and CE, the comparisons indicate that the effectiveness of the RK and merging methods on the higherelevation stations are similar to those of the RK and merging methods on the lower-elevation stations. For those stations with elevations above 1000 m but that are located outside or on the watershed boundary, the pair of averaged values for percentage reduction in RMSE and percentage improvement in CE from the RK and merging methods are 11.42% and 8.62% and 10.87% and 8.2%, respectively. Although the performances of the RK and merging approaches for the high-elevation stations outside or close to the watershed boundary are not as good as those for stations inside the watershed, their effectiveness in adjusting the radar rainfall data can still be considered to be reasonable.
To examine further the perfoance trend of the RK and merging methods, the overall averaged values for the percentage reduction in RMSE and percentage improvement in CE for the event of Typhoon Morakot are calculated. At stations that are inside the watershed boundary and with elevations over 1000 m, they are 40.93% and 38.5% and 38.81% and 37.28%, respectively, for the RK and merging methods. The pair values are 38.58% and 33.92% and 37.02% and 32.6% from RK and merging methods, respectively, at stations that are inside the watershed boundary but with elevations under 1000 m. Again, similar performance results can be found when applying both the RK and merging methods to either the higher-elevation or lower-elevation stations. For the 11 stations that are located outside (or on) the watershed boundary and with elevations over 1000 m, the averaged pair values of percentage reduction in RMSE and percentage improvement in CE are 15.6% and 22.2% and 11.22% and 16.47%, respectively, for the approaches of RK and merging. Again, the performances of RK and merging at those high-elevation stations are still acceptable.
With the above-described performance evaluations of the RK and merging methods on two major typhoon events, it is demonstrated that RK and merging are two approaches that can be effectively applied to correct the raw radar rainfall data reasonably well, even using the limited gauge data collected at higher-elevation stations. It is known that all the geospatial interpolation techniques rely heavily on the sampling quality. High-density gauge networks are challenging to set up in highelevation regions. Therefore, for practical application considerations, under inadequate sampling conditions, the RK and merging techniques would outperform other conventional methods due to their dual deterministic and stochastic processes.
Overall, the advantage of using multivariate interpolation techniques, such as RK and merging, is reflected by the demonstrated improvement made to the raw QPESUMS data. Both the RK and merging methods performed consistently in correcting data, especially at locations outside of the watershed boundary where sampling may be insufficient. With the detailed rainfall distribution information provided by the radar estimates of the QPESUMS, the RK and merging methods could be applied to update the rainfall reasonably well through the interpolation procedures. As a result, more realistic radar rainfall data can be generated for potential modeling uses.

Conclusions
Radar rainfall recorded for hydrological studies provides not only the time variation in an event but also the detailed spatial distribution over an area. While radar data preserve the spatial characteristics of the event, anomalies and uncertainty remain in the converted rainfall values. Before their use as inputs in modeling, proper adjustments to the radar rainfall values using gauge data are often required. In this study, the development and testing of two geostatistically based algorithms, regression-kriging (RK) and merging, for correcting radar rainfall data at the high-elevation Chenyulan river watershed in Taiwan were presented. Both of the techniques are multivariates through which the projected rainfall values were made by combining the gauge observations together with the correlated auxiliary variables-radar estimates in this case-as secondary predictors, such that the accuracy of the corrected rainfall data could be improved, especially at locations where the densities of the gauge distributions were inadequate.
On average, four typhoons with intensities of moderate or higher strike Taiwan every year. Severe property damage and the loss of lives are often inevitable. Reliable warning systems with accurate estimates of rainfall input are in great demand to provide emergency response agencies with vital information for early action planning in an effort to minimize any potential impact. Radar rainfall data, extracted from the QPESUMS system, are subject to correction using ground rain gauge observations.
The proposed algorithms, which were written in R, a statistical programming language, and model parameters were applied and verified with the events of five typhoons-namely typhoons Kalmaegi, Morakot, Fungwong, Sinlaku, and Fanapi-that hit Taiwan from 2008 to 2010. The interpolated rainfall values obtained with the RK and merging techniques were cross validated with the true gauge data and compared to the interpolated results from the ordinary kriging (OK) method, a univariate technique. The performance for each of the RK and merging methods was examined by analyzing the reduction percentage of the RMSE and the improvement percentage of the CE values from the comparisons between the corrected rainfall values and the original QPESUMS data. The plots of time variations in rainfall depth and data scattering along the 45-degree reference line were also generated to further the confirmation of the improvement made to the adjusted radar data. Based on the comparison results, it was clearly revealed that both the RK and merging algorithms could effectively produce reliable rainfall data for the study watershed, which is located at high elevations in a mountainous environment. At locations inside the watershed, both of the approaches were shown to be able to improve the rainfall values, with the RMSE reduced by up to nearly 50%. The improvement in the CE values of the corrected rainfall values was also evidenced by the data comparisons. Therefore, it was demonstrated that the original QPESUMS rainfall data could be adjusted with better accuracy for most stations tested. From the error analysis results, it was fair to indicate that the RK procedure, when applied to the five typhoon events, consistently made better adjustments on the QPESUMS rainfall data than the merging method did. In addition, the RK and merging methods were shown to outperform the OK method for correcting the radar data, especially for the locations that had an inadequate density of gauge station distribution. It is concluded that both the geostatistical based models, especially the RK method, can serve as practical and effective interpolation tools for reasonably correcting the raw radar rainfall data to values close to the gauge measurements. Since deep learning has been applied to disaster monitoring in recent years [65], for future work, employing recurrent neural networks and long short-term memory in geostatisticalbased models is worth trying, due to their ability to handle large amounts of time-sequence data such as radar data.