Comparison of spatial interpolation methods for the estimation of precipitation patterns at different time scales to improve the accuracy of discharge simulations

Precipitation is one of the most important variables in the water cycle. Accurate continuous spatial precipitation data plays a significant role in the hydrological planning, modeling, and management of water resource systems (Parker et al. 2019). As precipitation is often only recorded at individual gauge sites and is very dispersed at larger scales (Vicente-Serrano et al. 2003), continuous spatial data over a region cannot be directly collected from measurements.

In order to solve this problem, various spatial interpolation methods for the estimation of precipitation have been developed based on point measurements over the past several years (Chaubey et al. 1999; Zhang & Srinivasan 2009). Spatial interpolation methods range from conceptually simple weighting methods, such as Thiessen polygons (Thiessen 1911) or Inverse Distance Weighting schemes, to more complex approaches using artificial intelligence and covariates (Burrough & Mc Donnel 1998; Vicente-Serrano et al. 2003). In addition to elevation data, other parameters extracted from digital elevation models (DEMs) are used as covariates for interpolating precipitation. Grid data, such as Radar (NEXRAD) (Teegavarapu et al. 2012) and satellite data acquired by Tropical Rainfall Measuring Mission (TRMM) products (Wagner et al. 2012), are increasingly used in spatial interpolation as covariates. The mixed, or multi-method ensemble, approach has been increasingly used to improve model estimations (Hagedorn et al. 2005). The contribution of each individual model in the mixed, or multi-method ensemble, approach is represented by its weight. The use of simple weighting methods is still very common, as there are always insufficient or low-resolution covariate data and requires only low intensity computation. As both simple weighting methods and complex approaches have their own advantages and disadvantages based on specific case studies, it is necessary to select an appropriate spatial interpolation method for a given input dataset (Burrough & Mc Donnell 1998). Many comparative studies on precipitation interpolation techniques have been conducted (Vicente-Serrano et al. 2003; Li & Heap 2011; Bennet et al. 2013; Plouffe et al. 2015; Ding et al. 2018); however, there is little consensus on the optimal interpolator for precipitation (Dirks et al. 1998; Zimmerman et al. 1999; Price et al. 2000; Plouffe et al. 2015). Although some studies have compared the interpolation performance for different station densities (Goudenhoofdt & Delobbe 2009), very few have considered different time scales (Dirks et al. 1998; Bárdossy & Pegram 2014; Plouffe et al. 2015). As precipitation data properties differ with their time scales, the performance of a spatial interpolation method is affected. Therefore, it is important to develop a process for determining which method of spatial interpolation is best suited for different time scales.

The aims of this study are to provide a framework and useful suggestions on how to determine the most effective spatial interpolation methods for precipitation at different time scales, and to improve the accuracy of discharge simulations. The accuracy assessment for the precipitation interpolation results is done by applying cross-validation techniques (Hattermann et al. 2005) at different time scales, and then simulating the corresponding runoff driven by the precipitation interpolation results using a unifying hydrological model at the different time scales. The propagation of the uncertainties from the precipitation spatial interpolation method at each time scale is also discussed.

In order to evaluate the performance of various precipitation interpolation methods with regard to their suitability at different time scales, daily, monthly, and annual precipitation distributions were selected for testing via cross-validation schemes. The results from each precipitation interpolation method were also used to drive hydrological models in order to analyze the effects of the different precipitation interpolation methods on the discharge simulation. A framework of the entire workflow process employed in this study is depicted in Figure 1, with associated abbreviations shown in Table 1.

Spatial interpolation methods can be grouped into four categories: local methods, global methods, geostatistical methods, and mixed methods (Burrough & Mc Donnell 1998; Plouffe et al. 2015). Six different methods (Table 1) were applied and compared in our case study. These methods were evaluated on daily, monthly, and annual scales.

In order to eliminate the overreach of biasing inferences on or the bias from only one of several methods considered, simple linear regression is often used, and both the target and simulations are pre-processed using the Box-Cox transformation to align them more closely with the Gaussian distribution (Chen et al. 2015). The parameters θ and w can be estimated by the maximum the logarithm of the likelihood function shown in Equation (14). As the function cannot be maximized analytically, the Expectation-Maximization (EM) algorithm is used to numerically solve the problem (Wu 1983; McLachlan & Krishnan 1997; Raftery et al. 2005).

Wang & Tang (2014) and Zhao et al. (2016) derived a unifying catchment water balance model for different time scales through the maximum entropy production principle. The Budyko-type model at the annual scale, the ‘abcd’ model at the monthly scale, and the Soil Conservation Service (SCS) curve number model at the event scale are special cases of unifying catchment water balance models. In order to minimize the uncertainties of the hydrological model structures, the effects of the various interpolation inputs on the water balance and runoff dynamics at different time scales were assessed using the Soil & Water Assessment Tool (SWAT) based on the SCS model at the daily scale, ‘abcd’ at the monthly scale, and Budyko-type at the annual scale.

Since the SWAT is a continuous, semi-distributed, process-based basin scale model (Arnold et al. 1998), the basin was divided into multiple sub-basins which were then further subdivided into hydrologic response units (HRUs) (Wu et al. 2019). Water entering each sub-basin stream via direct runoff, soil lateral flow, or groundwater flow was then routed through the watershed stream network (Aliyari et al. 2019). Direct runoff was estimated using the Soil Conservation Service Curve Number (SCS-CN) model (Hjelmfelt 1991). Precipitation was partitioned into direct runoff and soil wetting, where soil wetting included initial abstraction and continuing abstraction. The direct runoff was assumed to equal the product of the remaining water after initial abstraction and the ratio of continuing abstraction to its potential continuing abstraction. Lateral flow, as a function of the interflow-recession coefficient, was calculated using the kinematic storage model (Neitsch et al. 2005; Chen et al. 2019). Groundwater flow in shallow (active) and deep (inactive) aquifers was simulated using a conceptual linear one-reservoir (shallow aquifer storage) approach (Fu et al. 2015; Pan et al. 2017).

The ‘abcd’ model is a nonlinear monthly hydrological model that was originally proposed by Thomas (1981) for national water assessment. The key assumption of the ‘abcd’ model is that the evapotranspiration opportunity is nonlinearly related to the available water in a manner such that the evapotranspiration opportunity increases quickly with the available water for water-limited conditions, but asymptotically approaches a maximum constant value for energy-limited conditions (Martinez & Gupta 2010). The structure of the model can be found in Thomas (1981), Alley (1984), Fernandez et al. (2000), and Sankarasubramanian & Vogel (2002).

At the watershed scale, if the water storage change in the mean annual water balance was negligible, the mean annual precipitation was portioned into runoff and evaporation (Budyko 1958). The portioning of the precipitation was determined by the competition between the available water (i.e., the annual precipitation) and available energy as measured by potential evaporation. The evaporative index, expressed by the ratio between the actual evaporation and precipitation, is a non-linear function of the climate dryness index, which was determined by the ratio of potential evaporation to precipitation. Then, the runoff was simply the difference between the precipitation and actual evaporation (Yang et al. 2008).

The study area is located upstream of Yichang city in Hubei and is defined as the upper reach of the Changjiang River (24.5 °N–27.8 °N; 90.6 °E–111.5 °E), which has a large topographic gradient ranging from 200 m to 6,500 m. The river is 4,500 km long with a basin area of 1.0 × 10⁶ km² (Figure 2). The upper reach of the Changjiang River has a subtropical, humid monsoon climate, which is characterized by seasonal rainfall from June to October. Owing to its extensive geographic area and various sources of precipitation (i.e., the South Asia summer monsoon (SASM), Plateau summer monsoon, and elevation), the precipitation spatial distribution is irregular at different time scales. As such, the precipitation spatial interpolations cannot ignore the impacts of time scales in the study area. The mean annual precipitation amount is 1,232 mm, and the annual average temperature is 15–18 °C. The dominant soil types are purple, yellow, and paddy soils. The Three Gorges Dam, the world's largest hydropower project to date, is located 44 km upstream of the Yichang hydrological gauge station.

Daily precipitation measurements at 53 gauge stations, with continuous data collected from 1,960 to 2008, were collected from the website of the China Meteorological Bureau (http://data.cma.cn/). These precipitation stations displayed a relatively uniform distribution over the study area as shown in Figure 2. The daily discharge data at the Yichang gauge station from 1960 to 2008 were provided by the Bureau of Hydrology of the Changjiang Water Resources Commission of China.

The critical problem associated with the TSA is selection of the optimal polynomial order (Chayes 1970). In this study, the degree of fitting (AIC) was used to determine the polynomial order of the TSA. The lower the AIC value a particular method produced, the better fitness the observed data provided. As presented in Table 2, a third order polynomial was assumed to be sufficient to capture all the daily, monthly, and annual precipitation variations.

To identify the optimal parameters of the IDW method at the daily, monthly, and annual time scales, the α value, number of sample points, N, the research radius, the minimum MAE and RMSE, and the maximum R should be determined. The values of MAE, RMSE, and R, shown in Figure 3, revealed that the optimal values for parameters N and α can be determined as shown in Table 3.

Three semivariogram (Gaussian, spherical, and exponential) models were adopted for fitting the empirical model in the OK and CK spatial interpolation methods. The cross-validation technique was used to calibrate the parameters of the semivariogram models. The best-fitting method was determined by selecting the model with the lowest MAE and RMSE, or the highest R. There were no significant differences among the three semivariogram models or between the OK and CK methods at the daily (Figure 4(a)), monthly (Figure 4(b)), and annual (Figure 4(c)) time scales. However, the spherical semivariogram model was the best, as determined by the average statistical criteria values (shown in Table 4) and was selected by the OK and CK spatial interpolation methods.

The global optimal BMA weights for the log-likelihood function (13) were obtained by using the EM algorithm with daily, monthly, and annual time steps (shown in Figure 5). As there were five models, the average weights were approximately 0.2. The uncertainty of the daily weights was the highest (Figure 5(a)), while the annual was the lowest (Figure 5(c)). Only the expected BMA estimation was adopted in this study.

After the parameters of the five interpolation methods (TSA, OK, CK, TSA-OK, and BMA) were calibrated, their performance, including NN and IDW, without calibration parameters was evaluated using the goodness-of-fit indicators (MAE, R, and RMSE) based on the cross-validation scheme. Considering the observed data from the entire basin, the performance of all seven spatial interpolation methods is shown in boxplots in Figure 6. The mixed methods (TSA-OK and BMA) performed best at the daily, monthly, and annual time scales. BMA encompassed five individual methods, while TSA-OK encompassed only two single methods. Their performance was similar considering their MAE, R, and RMSE values at the monthly and annual time scales. The BMA weights at the daily scale (shown in Figure 5) were more dispersive than at the monthly and annual scales. It can be seen that the performance of BMA combined with seven methods was better than TSA-OK at the daily scale. The more dispersive the BMA weights, the more accurate the results from the methods that combined one or more approaches. If the weights of every single method in BMA were equal or similar, there was no significant difference between the mixed methods that were combined with various numbers of single methods. As the spatial precipitation correlations for short periods (daily) varied more than those for longer periods (monthly and annually), it was determined that a mixed method with more single models was more accurate than a mixed method with fewer ones. Although the best performance of the mixed method might be partly due to the different number of parameters, the input for every method was just the precipitation series, and recommendation of a mixed method will still be useful in practice.

Of the single methods evaluated, the best performance was from IDW. OK was better than CK, TSA was better than NN and CK, and OK was better than TSA at the monthly and annual time scales, while the former was worse than the latter at the daily scale. The spatial variation of the monthly precipitation was higher than that of the annual precipitation. The values of R increased from the daily, monthly, and annual time scales, while their variation ranges decreased, which was different from the values of R at each site. The performance of each single method was not significantly different at the monthly and annual scales from the perspective of the entire study area. The annual data compensated for the possible differences in daily precipitation amounts, as in the case of moving storms, which was one reason to improve the interpolation methods. Therefore, while any method could be recommended to estimate the monthly or annual spatial precipitation distribution, IDW was the primary recommendation.

The performance of all seven methods was also examined at 53 gauge stations. At a given station, the larger MAE, R, and RMSE values are shown in darker red, while the smaller values are shown in darker blue, as seen in Figure 7. According to the gauge station distribution shown in Figure 2, the better performance locations were found downstream of the study area where the elevation was lower. Estimations for the upstream sites from the 3rd to the 30th gauge stations were not as good as those downstream. As the precipitation gauge station density was low at the border of the study area (see Figure 3), the poor performance at these sites can be partly attributed to the lack of precipitation gauges with similar characteristics. Considering the R values, the performance of each method at the monthly and annual time scales improved compared to the daily time scale (e.g., No. 13 gauge site). However, the R values at the monthly scale were better than those at the annual scale for most sites. As precipitation in the study area comes mainly from monsoons, the spatial correlations between sampling sites were strong and beneficial to spatial interpolation. The performance of each method at some sites was poor for every time scale (e.g., No. 17 gauge site).

In order to directly compare the performance of each method at a given gauge site, the worst method was given a value of ‘one’ according to the MAE, R, and RMSE values at a particular time scale, while the best method was given a value of ‘seven’. Applying these ‘points’ to each gauge site, Figure 7 is transformed into Figure 8. It can be easily seen that the performance of each method was different for different time scales. The mixed methods (TSA-OK and BMA) performed best at the daily scale, while CK, OK, and NN performed worst among the methods. As the time scales increased to the monthly and annual scales, no method with significant strengths emerged. However, some models demonstrated good performance at some sites, but poor performance at other sites. For example, BMA for monthly precipitation produced the lowest MAE and RMSE values at the second gauge site, but produced higher values at the 17th site. The rank score color as determined by the MAE, R, and RMSE values were not always the same at a given gauge site for a particular time scale. It should be noted that the performance of the NN method was high at the monthly and annual time scales, while it demonstrated poor performance at the daily time scale.

The performance of the SWAT, ‘abcd’, and Budyko hydrological models driven by the interpolated precipitation were evaluated using the observed discharge data at the Yichang gauge station from 1960 to 2008 at the daily, monthly, and annual time scales. The ranks of the NSE values for the discharge simulation at the daily, monthly, and annual time scales were the same as that of the R values for the spatial interpolation methods at the 53 gauge stations (shown in Table 5 and Figure 7). However, the ranks of PBIAS were consistent with those of the R values with respect to the entire study area (shown in Table 5 and Figure 6). The highest NSE values were obtained from the ‘abcd’ model at the monthly scale for the calibration and validation periods. The lowest PBIAS values were obtained from the Budyko model at the annual scale for the calibration and validation periods. Since the distribution of precipitation at the monthly scale impacts the discharge at the Yichang gauge station, the annual discharge was primarily determined by the accuracy of the interpolated precipitation over the entire study area. A more favorable NSE value indicated that the BMA spatial interpolation method provided better results except for the modeled runoff at the daily time scale for the calibration period. Additionally, the ranks of the performance of the discharge simulation at the monthly and annual time scales were consistent with the ranks of precipitation estimation with respect to the NSE values for the calibration and validation periods. The differences in the discharge simulation at the daily scale, driven by the spatially interpolated precipitation, were significant with respect to the NSE or PBIAS values (shown in Table 5). The range of their differences was 0.05–0.06 for NSE (which accounts for nearly 10% of the average NSE value) and 5.50–5.76% for PBIAS (which accounts for nearly 50% of the average PBIAS value). However, the differences at the monthly and annual time scales were smaller. If the metrics NSE and PBIAS were determined for the total period from 1960 to 2008, BMA would be found to have the best performance for the discharge simulation at the three time scales. However, the ranks of the discharge performance driven by another spatial precipitation interpolation scheme were different from those of spatial precipitation based on cross-validation techniques. It can also be inferred from the discharge simulation shown in Table 5 that the sample splitting for calibration and validation led to uncertainties in the discharge simulation.

Scatterplots depicting the observed discharge versus the simulated discharge at different time scales for low discharges exhibited a closer linear relationship than at high discharges. In addition, the simulated discharge was lower than the observed at the daily and monthly time scales for the calibration and validation periods. The observed vs. simulated discharge plots driven by the NN spatially interpolated precipitation at the daily scale for the validation period in Figure 9(a₂) showed more concordance between the simulated and observed values than in the calibration period (Figure 9(a₁)). The spatial precipitation relationships between each gauge station at the daily scale were not as strong as at the monthly and annual time scales. The land surface conditions, which consider such things as water storage in the soil, contributed to the flow generation and routing. Therefore, the performance of the discharge simulation at the daily scale was not exactly the same as the corresponding precipitation estimate. The difference in large discharge values for the calibration period was more than in the validation period, while the NSE and PBIAS values for the calibration period were less than in the validation period (Table 5).

Some of the error associated with the results of the areal precipitation estimation and discharge simulation can be accounted for when looking at the entire study area. Given that the upstream portion of the study area had significantly fewer available sampling points than the downstream portion (Figure 2), the precipitation discrepancies become problematic, as precipitation for large regions that cover many different elevations was estimated with little input data. Finding accurate spatial relationships for precipitation between sparse and relatively dense gauge stations for different time scales that could improve spatial precipitation estimation and discharge simulation was the aim of this study.

It can be concluded that mixed methods, such as TSA-OK and BMA, performed best at the daily, monthly, and annual time scales. If the weights of every single method in BMA were equal or similar, there was not much difference between the mixed methods, which were combined with various single methods. As a single method, IDW produced low values for the MAE and RMSE metrics for all three time scales. However, it did not perform the same at all sites. As the time scale increased to monthly and annual, there were no methods with significant strengths. However, several models showed good performance at some sites, but poor performance at other sites. There was no method that performed well at every site in all three time scales. As the time scale increased, the differences in interpolated precipitation between each method shrank. The selection of spatial interpolation methods for areal precipitation at small time scales should be given extra consideration. The mixed methods demonstrated a significant advantage in spatially interpolating precipitation for small time scales.

The adopted hydrological models for the discharge simulation should be different from their time scales. SWAT based on SCS, ‘abcd’, and the Budyko hydrological models have been proven as unifying catchment water balance models for different time scales through the maximum entropy production principle (Wang & Tang 2014; Zhao et al. 2016). Propagation of the uncertainty of spatial precipitation interpolation methods were comparable at different time scales. It was found that the uncertainty caused by the different spatial precipitation interpolation methods decreased with the increase in the time scale with respect to the discharge simulation. The performance of the daily discharge simulation was not consistent with the performance of the daily precipitation interpolation, because factors such as water storage and initial water conditions influence the discharge simulation in addition to precipitation. It was typical for the discharge simulation to underestimate high discharge. Ultimately, there was no dominant spatial precipitation interpolation for the discharge simulation at small time scales (i.e., daily). Although the selection of spatial interpolation methods for areal precipitation estimation and discharge simulation at different scales have been proposed based on our case study, more case studies should be done for more general conclusion and guideline.

The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (Nos. 51579183, 91647106, 51879194, and 51525902). This work is also partly funded by the Ministry of Foreign Affairs of Denmark and administered by Danida Fellowship Centre (File number: 18-M01-DTU).