Abstract

Interpolating precipitation data is of prime importance to hydrological design, modeling, and water resource management. Various models have been developed that estimate spatial precipitation patterns. The purpose of this study is to analyze different precipitation interpolation schemes at different time scales in order to improve the accuracy of discharge simulations. The study was carried out in the upstream area of the Changjiang River basin. The performance of all selected methods was assessed using cross-validation schemes, with the mixed methods ultimately displaying the best performance at all three time scales. However, the differences in performance between the spatial interpolation methods decreased with increasing time scales. The unifying catchment Soil and Water Assessment Tool (SWAT), ‘abcd’, and the Budyko equation were employed at the daily, monthly, and annual scales, respectively, to simulate discharge. The performance of the discharge simulation at the monthly and annual time scales was consistent with their ranks of spatial precipitation estimation. For coarse, or long period, precipitation, there were no significant differences. However, the mixed methods performed better than the single model for the daily, or short, time scale with respect to the accuracy of the discharge simulation.

HIGHLIGHTS

  • The performances of different precipitation interpolation schemes at different time scales have been established.

  • The mixed spatial interpolation schemes for precipitation are found to be the best at all time scales due to its own parameters and structure.

  • The accuracy of discharge simulation is found to be sensitive to the spatial interpolation schemes for precipitation at short time scale while there are no significant differences at monthly and annual time scales.

  • The importance of time scales in the precipitation interpolation are recommended to be considered before determining the areal precipitation.

  • The proposed framework for selecting the methods can also be applied to other spatial precipitation interpolation methods.

INTRODUCTION

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.

METHODOLOGY

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.

Table 1

Spatial interpolation methods for estimating precipitation distribution

Spatial interpolation methodsGroupAbbreviation
Nearest Neighbor Local NN 
Inverse Distance Weighting Local IDW 
Trend Surfaces and Regression Global TSA 
Ordinary Kriging Geostatistical OK 
Co-kriging Geostatistical CK 
TSA combined with OK Mixed TSA-OK 
Bayesian Model Averaging Mixed BMA 
Spatial interpolation methodsGroupAbbreviation
Nearest Neighbor Local NN 
Inverse Distance Weighting Local IDW 
Trend Surfaces and Regression Global TSA 
Ordinary Kriging Geostatistical OK 
Co-kriging Geostatistical CK 
TSA combined with OK Mixed TSA-OK 
Bayesian Model Averaging Mixed BMA 
Figure 1

Framework illustrating steps taken to compare the accuracies of precipitation spatial interpolation methods and their effects of hydrological processes at different time scales.

Figure 1

Framework illustrating steps taken to compare the accuracies of precipitation spatial interpolation methods and their effects of hydrological processes at different time scales.

Spatial interpolation methods for estimating precipitation distribution

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.

Nearest Neighbor (NN) and Inverse Distance Weighting (IDW) are local spatial interpolation methods that utilize the concept of Tobler's first law to estimate precipitation by distance (Chen & Liu 2012). The NN (which is also called the Thiessen polygons interpolation method) method assumes that point pu at the unsampled spatial location is equal to the closest available rainfall, pi. If λu,i is the distance between locations u and i, the weight of λu,i is 1 for the location i corresponding to the minimum distance, du,i, and zero for all other locations (Thiessen 1911; Nikolopoulos et al. 2015). The IDW method, recommended by the Handbook of Hydrology by ASCE (1996), is easy to implement (Teegavarapu & Elshorbagy 2005). IDW (Shepard 1968) assumes that weights, λu,i, are also related to the distance, du,I, from the sampled point, pi, to the estimated point, pu, as given by Equation (1): 
formula
(1)
where N is the number of sample points, and α is the power parameter in the interpolation. The variable α is the friction distance (Vieux 2001), ranging from 1.0 to 6.0, and is the main factor affecting the estimation. It also controls the decay of the weights with distance. The higher the value of α, the greater the influence of the points closest to the interpolated location (Nikolopoulos et al. 2015).
Different from the basis of Tobler's first law, where the precipitation is estimated by distance for the two local spatial interpolation methods (NN and IDW), trend surface analysis (TSA) is a global method that treats each value in a geographic area as a centroid in a grid pattern (Wang 2006). If U and V are geographic coordinates, pi,j is a single observation of the mapped variable at coordinates (Ui, Vj). The pu at an unmeasured spatial location can be estimated by Equation (2) (Wong et al. 2003): 
formula
(2)
where λi,j is the polynomial coefficient and N is the degree of the polynomial that represents the best surface of the precipitation pattern. λi,j and N are determined by using the least-squares approximation (Krumbein 1959) or minimum Akaike Information Criterion (AIC) (Bumhan & Anderson 2002; Soares et al. 2015). The major issue with TSA is selection of the appropriate functional form to model the trend (Sullivan & Unwin 2003).
As TSA, NN, and IDW do not exploit the statistical properties of the observational samples (Foehn et al. 2018), the weights of the Ordinary Kriging (OK) and Co-kriging (CK) methods are optimized based on information that is inherent in the measured data (Wagner et al. 2012). The weights in OK are obtained by solving Equations (3)–(5): 
formula
(3)
where γ(xixj) represents the value of the semivariogram function for the distance between points xi and xj. γ(xjxu) is the value of the distance between xj and the estimated location xu, and ϕ is the Lagrange parameter, which is equal to the mean distance. On the basis of the second-order stationary hypothesis, the semivariogram is derived by fitting a semivariogram model to the empirical one, and can be determined for all distances, h, by solving Equation (4): 
formula
(4)
Gaussian (Equation (5a)), spherical (Equation (5b)), and exponential (Equation (5c)) semivariogram models are employed to fit the empirical model shown in Equation (4): 
formula
(5a)
 
formula
(5b)
 
formula
(5c)
where C0, C0 + C, and a are the parameters known as the nugget, the sill, and the range, respectively. CK offers the possibility of considering more than one variable in the kriging interpolation, especially elevation (Pardo-Iguzquiza et al. 2011), and is expressed as follows: 
formula
(6)
where the subscripts ‘1’ and ‘2’ refer to the primary and secondary networks, respectively, and N1 and N2 are the number of available rain gauges in the primary and secondary networks over the period, respectively. Weights λ1 and λ2 are obtained by solving Equation (7). 
formula
(7)
In order to combine the physical relationships between climatic data and geographic or topographic variables, as done in the global TSA method, and the spatial correlation between the information recorded at weather stations, as done in OK, a mixed method, called TSA-OK, has been proposed to improve the accuracy of spatial interpolation methods. The TSA is an inexact interpolator, and has residual errors (as shown in Equation (8)). The residuals from the TSA method are then interpolated by the OK method (Equation (9)). The exact climatic variable for the final prediction is obtained through Equation (10): 
formula
(8)
 
formula
(9)
 
formula
(10)
Although all four categories of spatial interpolation methods have been widely used in the estimation of precipitation distribution, none has been able to select the correct method with certainty a priori. The uncertainties associated with selection of the method should be taken into account in order to reduce prediction error. Bayesian Model Averaging (BMA) addresses this issue by its inherent incorporation of method uncertainty into predictors via considering all competing methods (including variations in their parameters), and making predictions using a weighted average of all of them, with weights representing the posterior model probabilities (Schneider & Yaşar 2016; Zhang & Taflanidis 2019). Consider a predicted or interpolated precipitation, p, the corresponding target, or training data, pT, and an ensemble of all spatial interpolation method (i.e., TSA, NN, IDW, OK, and CK) simulations {f1, f2, f3, f4, f5} of the precipitation. According to the law of total probability (Bulmer 1979), the probability density function of p can be expressed as: 
formula
(11)
where Pr(p|fk) is the probability prediction given by the simulation of fk alone, and Pr(fk|pT) is the posterior probability of the method prediction, given the data set, pT. Pr(fk|pT) is also identified as a fractional statistical weight, wk, of the kth spatial interpolation method, whose magnitude reflects how well fk matches pT, and . In order to determine the prior Pr(p|fk) for the parameters and the prior probability of each method, it is computationally convenient to assume that Pr(p|fk) is a Gaussian distribution defined by a mean, μk, and variance, σk2 (Raftery et al. 2005; Duan & Phillips 2010). 
formula
(12)
where g(•) is the Gaussian probability density function, and θk = [μk, σk2] is the parameter vector.
Substituting (12) into (11), 
formula
(13)
The log-likelihood function, l, was maximized to obtain both the Bayesian weight, wk, and the parameter vector, θk, for its convenient computation. The log-likelihood function, l, is approximated as: 
formula
(14)

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).

Hydrological models

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.

SWAT model

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).

 ‘abcd’ model

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).

Budyko model

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).

Statistical criteria for assessing the agreements of the spatial interpolation methods and hydrological models

Statistical criteria for the spatial interpolation used to estimate precipitation

The performance of every spatial interpolation used to estimate precipitation was assessed and compared using the leave-one-out, or cross-validation, approach (Goovaerts 2000). This approach estimated precipitation at a gauge location using all observations except the one corresponding to the interpolation location (Isaaks & Srivastava 1989; Foehn et al. 2018). The comparison criteria were the Mean Absolute Error (MAE), Pearson Product moment correlation measures (R), and Root-Mean-Square Error (RMSE), which measure overall accuracy (MAE and R) or reflect the extreme outliers (RMSE) (Berndt & Haberlandt 2018). The procedure was undertaken on different time scales (daily, monthly, and annually in this study) for each of the spatial interpolation methods. 
formula
(15)
 
formula
(16)
 
formula
(17)
where Oi and are the observation and estimation at the ith gauge station, respectively.

Performance assessment of hydrological models

To evaluate the fitness of the hydrological models, Nash-Sutcliffe Efficiency (NSE) (Nash & Sutcliffe 1970) and Percentage Bias (PBIAS) (Moriasi et al. 2007) were used to measure the simulation accuracy of the discharge and are given by Equations (18) and (19): 
formula
(18)
 
formula
(19)

STUDY AREA AND DATA SET

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 × 106 km2 (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.

Figure 2

Location of precipitation gauges used for assessment of interpolation methods in upstream of Changjiang River basin.

Figure 2

Location of precipitation gauges used for assessment of interpolation methods in upstream of Changjiang River basin.

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.

RESULTS AND DISCUSSION

Parameters of precipitation spatial interpolation calibration

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.

Table 2

The AIC values for the first four orders of TSA

N*DailyMonthlyAnnual
N = 1 73.63 358.08 596.35 
N = 2 70.55 349.33 589.44 
N= 3 69.57 348.51 587.81 
N = 4 70.41 350.2 588.95 
N*DailyMonthlyAnnual
N = 1 73.63 358.08 596.35 
N = 2 70.55 349.33 589.44 
N= 3 69.57 348.51 587.81 
N = 4 70.41 350.2 588.95 

*N is the polynomial degree (order) of TSA.

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.

Table 3

The determined parameters for the IDW spatial interpolation methods

ParameterDailyMonthlyAnnual
N 
α 1.6 1.8 1.6 
ParameterDailyMonthlyAnnual
N 
α 1.6 1.8 1.6 
Figure 3

The values of MAE, R and RMSE with the variations of the numbers of rainfall stations N and the α values at (a) daily, (b) monthly and (c) annual time scales.

Figure 3

The values of MAE, R and RMSE with the variations of the numbers of rainfall stations N and the α values at (a) daily, (b) monthly and (c) annual time scales.

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.

Table 4

The average statistical criteria values of the semivariogram models for OK and CK spatial interpolation methods

Semivariogram modelsStatistical criteriaOK
CK
DailyMonthlyAnnualDailyMonthlyAnnual
Gaussian MAE (mm) 2.76 27.32 164.55 2.95 28.33 173.39 
R 0.44 0.67 0.77 0.40 0.66 0.70 
RMSE (mm) 4.26 37.37 212.56 4.48 38.38 242.65 
Spherical MAE (mm) 2.74 27.28 164.48 2.90 28.12 173.27 
R 0.45 0.68 0.78 0.42 0.67 0.72 
RMSE (mm) 4.25 37.34 212.24 4.41 38.32 242.55 
Exponential MAE (mm) 2.77 27.31 165.48 2.97 28.23 173.31 
R 0.45 0.67 0.75 0.41 0.67 0.71 
RMSE (mm) 4.27 37.39 212.59 4.46 38.40 243.08 
Semivariogram modelsStatistical criteriaOK
CK
DailyMonthlyAnnualDailyMonthlyAnnual
Gaussian MAE (mm) 2.76 27.32 164.55 2.95 28.33 173.39 
R 0.44 0.67 0.77 0.40 0.66 0.70 
RMSE (mm) 4.26 37.37 212.56 4.48 38.38 242.65 
Spherical MAE (mm) 2.74 27.28 164.48 2.90 28.12 173.27 
R 0.45 0.68 0.78 0.42 0.67 0.72 
RMSE (mm) 4.25 37.34 212.24 4.41 38.32 242.55 
Exponential MAE (mm) 2.77 27.31 165.48 2.97 28.23 173.31 
R 0.45 0.67 0.75 0.41 0.67 0.71 
RMSE (mm) 4.27 37.39 212.59 4.46 38.40 243.08 
Figure 4

The performance boxplots of the three semivariogram models for OK and CK methods at 53 precipitation gauge sites: (a) daily, (b) monthly and (c) annual. The whiskers of the boxplots are drawn to the 90th and 10th percentiles. The center line (red line in the rectangle) is the median and the dot is the average values. The rectangle is defined by the 75th and 25th percentiles as upper and lower hinges (same meanings in the following boxplots Figures 5 and 6). Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

Figure 4

The performance boxplots of the three semivariogram models for OK and CK methods at 53 precipitation gauge sites: (a) daily, (b) monthly and (c) annual. The whiskers of the boxplots are drawn to the 90th and 10th percentiles. The center line (red line in the rectangle) is the median and the dot is the average values. The rectangle is defined by the 75th and 25th percentiles as upper and lower hinges (same meanings in the following boxplots Figures 5 and 6). Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

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.

Figure 5

BMA weights boxplots reflecting the matching degrees for NN, IDW, ISA, OK, CK spatial interpolation methods: (a) daily; (b) monthly; (c) annual.

Figure 5

BMA weights boxplots reflecting the matching degrees for NN, IDW, ISA, OK, CK spatial interpolation methods: (a) daily; (b) monthly; (c) annual.

Performance of the spatial interpolation methods for estimating precipitation at different time scales

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.

Figure 6

The performances boxplots of selected seven precipitation interpolation methods for study area at (a) daily; (b) monthly; (c) annual scales.

Figure 6

The performances boxplots of selected seven precipitation interpolation methods for study area at (a) daily; (b) monthly; (c) annual scales.

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).

Figure 7

The performances of spatial interpolation methods at 53 gauge stations according to MAE, R and RMSE at daily (a), monthly (b) and annual (c) time scales. The site ID indicates 53 precipitation gauge sites (see Figure 2). The color bar indicates the values of metrics, and the blue and red represent the worst and best model performance, respectively. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

Figure 7

The performances of spatial interpolation methods at 53 gauge stations according to MAE, R and RMSE at daily (a), monthly (b) and annual (c) time scales. The site ID indicates 53 precipitation gauge sites (see Figure 2). The color bar indicates the values of metrics, and the blue and red represent the worst and best model performance, respectively. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

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.

Figure 8

The performance rank of seven spatial interpolation methods at 53 gauge stations according to MAE, R and RMSE at daily, monthly and annual time scales. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

Figure 8

The performance rank of seven spatial interpolation methods at 53 gauge stations according to MAE, R and RMSE at daily, monthly and annual time scales. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/nh.2020.146.

Performance of hydrological models driven by precipitation at different time scales

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.

Table 5

Performances of hydrological models driven by differently interpolated precipitation inputs

SWAT (Daily)
abcd (Monthly)
Budyko (Annual)
NSE
PBIAS(%)
NSE
PBIAS(%)
NSE
PBIAS(%)
Spatial precipitation interpolation schemeIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
NN 0.702 0.744 0.699 − 13.05 −8.91 −13.54 0.916 0.908 0.914 7.95 4.72 6.99 0.541 0.824 0.692 −0.15 −0.58 −0.40 
IDW 0.749 0.729 0.828 − 8.74 −11.36 −9.57 0.936 0.916 0.930 − 1.17 −5.34 −2.39 0.602 0.838 0.728 0.09 −0.45 −0.16 
TSA 0.732 0.723 0.728 − 9.98 −11.82 −10.59 0.933 0.909 0.926 − 2.02 −6.57 −3.33 0.569 0.825 0.706 0.10 −0.45 −0.15 
OK 0.724 0.718 0.723 − 10.96 −11.91 −11.15 0.930 0.908 0.923 − 2.15 −5.46 −3.08 0.560 0.822 0.700 0.12 −0.51 −0.17 
CK 0.701 0.721 0.696 − 12.29 −11.58 −13.04 0.920 0.906 0.916 0.27 −2.03 −0.35 0.551 0.821 0.695 0.22 −0.51 −0.10 
TSA-OK 0.736 0.692 0.732 − 9.90 −8.91 −10.51 0.936 0.908 0.930 1.53 4.72 −.74 0.633 0.824 0.750 0.14 0.49 −0.14 
BMA 0.753 0.684 0.831 − 7.55 11.36 7.98 0.939 0.916 0.932 − 1.15 −5.34 −2.17 0.649 0.838 0.760 0.09 −0.37 −0.12 
SWAT (Daily)
abcd (Monthly)
Budyko (Annual)
NSE
PBIAS(%)
NSE
PBIAS(%)
NSE
PBIAS(%)
Spatial precipitation interpolation schemeIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
NN 0.702 0.744 0.699 − 13.05 −8.91 −13.54 0.916 0.908 0.914 7.95 4.72 6.99 0.541 0.824 0.692 −0.15 −0.58 −0.40 
IDW 0.749 0.729 0.828 − 8.74 −11.36 −9.57 0.936 0.916 0.930 − 1.17 −5.34 −2.39 0.602 0.838 0.728 0.09 −0.45 −0.16 
TSA 0.732 0.723 0.728 − 9.98 −11.82 −10.59 0.933 0.909 0.926 − 2.02 −6.57 −3.33 0.569 0.825 0.706 0.10 −0.45 −0.15 
OK 0.724 0.718 0.723 − 10.96 −11.91 −11.15 0.930 0.908 0.923 − 2.15 −5.46 −3.08 0.560 0.822 0.700 0.12 −0.51 −0.17 
CK 0.701 0.721 0.696 − 12.29 −11.58 −13.04 0.920 0.906 0.916 0.27 −2.03 −0.35 0.551 0.821 0.695 0.22 −0.51 −0.10 
TSA-OK 0.736 0.692 0.732 − 9.90 −8.91 −10.51 0.936 0.908 0.930 1.53 4.72 −.74 0.633 0.824 0.750 0.14 0.49 −0.14 
BMA 0.753 0.684 0.831 − 7.55 11.36 7.98 0.939 0.916 0.932 − 1.15 −5.34 −2.17 0.649 0.838 0.760 0.09 −0.37 −0.12 

⃰’I’ is for the calibration period (1960–1992); ‘II’ is for validation period (1993–2008) and ‘II’ is for whole period (1960–2008).

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(a2) showed more concordance between the simulated and observed values than in the calibration period (Figure 9(a1)). 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).

Figure 9

Scatterplots of the observed discharge versus the simulated mean discharge driven by different spatial interpolated precipitation at (a): daily, (b) monthly, (c) annual time scales. Subscript ‘1’ is from the calibration period while subscript ‘2’ is from the validation period. The solid line is the 1:1 line.

Figure 9

Scatterplots of the observed discharge versus the simulated mean discharge driven by different spatial interpolated precipitation at (a): daily, (b) monthly, (c) annual time scales. Subscript ‘1’ is from the calibration period while subscript ‘2’ is from the validation period. The solid line is the 1:1 line.

CONCLUSION

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.

ACKNOWLEDGEMENT

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).

REFERENCES

REFERENCES
Aliyari
F.
Bailey
R. T.
Tasdighi
A.
Dozier
A.
Arabi
M.
Zeiler
K.
2019
Coupled SWAT-MODFLOW model for large-scale mixed agro-urban river basins
.
Environ. Modell. Software
115
,
200
210
.
doi:10.1016/j.envsoft.2019.02.014
.
Arnold
J. G.
Srinivasan
R.
Muttiah
R. S.
Williams
J. R.
1998
Large area hydrologic modeling and assessment Part I: model development
.
J. Am. Water Resour. Assoc.
34
(
1
),
73
89
.
doi:10.1111/j.1752-1688.1998.tb05961
.
ASCE
1996
Hydrology of Handbook, Task Committee on Handbook of Hydrology (Management Group D of ASCE)
, 2nd edn.
American Society of Civil Engineers, ASCE
,
NewYork
.
Bennet
N. D.
Croke
B. F. W.
Guariso
G.
Guillaume
J. H. A.
Hamilton
S. H.
Jakeman
A. J.
Marsili-Libelli
S.
Newham
L. T. H.
Norton
J. P.
Perrin
C.
Pierce
S. A.
Robson
B.
Seppelt
R.
Voinov
A. A.
Fath
B. D.
Andreassian
V.
2013
Characterising performance of environmental models
.
Environ. Modell. Software
40
,
1
20
.
Berndt
C.
Haberlandt
U.
2018
Spatial interpolation of climate variables in Northern Germany – influence of temporal resolution and network density
.
J. Hydrol.: Reg. Stud.
15
,
184
202
.
Blumer
M. G.
1979
Principles of Statistics
.
Dover
,
New York
, p.
252
.
Budyko
M. I.
1958
The Heat Balance of the earth's surface
.
translated from Russian by N.A. Stepanova
.
U.S. Dep. of Commer.
,
Washington, DC
, p.
259
.
Bumham
K. P.
Anderson
D. R.
2002
Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach
.
Springer
,
Heidelber
.
Burrough
P.
Mc Donnell
R.
1998
Principles of Geographical Information Systems
.
Oxford University Press
,
New York
.
Chaubey
I.
Haan
C. T.
Salisbury
J. M.
Grunwald
S.
1999
Quantifying model output uncertainty due to the spatial variability of rainfall
.
J. Am. Water Resour. Assoc.
35
(
5
),
1113
1123
.
Chen
Y.
Yuan
W. P.
Xia
J. Z.
Fisher
J. B.
Dong
W. J.
Zhang
X. T.
Liang
S. L.
Ye
A. Z.
Cai
W. W.
Feng
J. M.
2015
Using Bayesian model averaging to estimate terrestrial evapotranspiration in China
.
J Hydrol.
528
,
537
549
.
Duan
Q.
Phillips
T. J.
2010
Bayesian estimation of local signal and noise in multimodel simulations of climate change
.
J. Geophys. Res.
115
(
D18
),
D18123
.
http://dx.doi.org/10.1029/2009JD013654
.
Fernandez
W.
Vogel
R. M.
Sankarasubramanian
A.
2000
Regional calibration of a watershed model
.
Hydrol. Sci. J.
45
(
5
),
689
707
.
Goudenhoofdt
E.
Delobbe
L.
2009
Evaluation of radar-gauge merging methods for quantitative precipitation estimates
.
Hydrol. Earth Syst. Sci.
13
,
195
203
.
Hagedorn
R.
Doblas-reyes
F. J.
Palmer
T. N.
2005
The rational behind the success of multi-model ensembles in seasonal forecasting-I. Basic concept
.
Tellus
57A
,
219
233
.
Hjelmfelt
A. T.
1991
Investigation of curve number procedure
.
J. Hydraul. Eng.
17
(
6
),
725
735
.
Isaaks
E. H.
Srivastava
R. M.
1989
An Introduction to Applied Geostatistics
.
Oxford University Press
,
New York
, pp.
351
368
McLachlan
G. J.
Krishnan
T.
1997
The EM Algorithm and Extensions
.
Wiley
,
New York
.
Moriasi
D. N.
Arnold
J. G.
Van Liew
M. W.
Bingner
R. L.
Harmel
R. D.
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. ASABE
50
(
3
),
885
900
.
Neitsch
S. L.
Arnold
J. G.
Kiniry
J. R.
Williams
J. R.
2005
Soil and Water Assessment Tool Theoretical Documentation, Version 2005
.
Pan
S. H.
Liu
D. D.
Wang
Z. L.
Zhao
Q.
Zou
H.
Hou
Y. K.
Liu
P.
Xiong
L. H.
2017
Runoff responses to climate and land use/cover changes under future scenarios
.
Water
9
(
7
),
475
.
https://doi.org/10.3390/w9070475
.
Pardo-Iguzquiza
E.
Rodríguez-Galiano
V. F.
Chica-Olmo
M.
Atkinson
P. M.
2011
Image fusion by spatially adaptive filtering using downscaling cokriging
.
ISPRS. J. Photogramm.
66
,
337
346
.
Parker
S. R.
Adams
S. K.
Lammers
R. W.
Stein
E. D.
Bledsoe
B. P.
2019
Targeted hydrologic model calibration to improve prediction of ecologically-relevant flow metrics
.
J. Hydrol.
573
,
546
556
.
Price
D. T.
McKenney
D. W.
Nalder
I. A.
Hutchinson
M. F.
Kesteven
J. L.
2000
A comparison of two statistical methods for spatial interpolation of Canadian monthly mean climate data
.
Agric. Meteorol.
101
(
2–3
),
81
94
.
Raftery
A. E.
Gneiting
T.
Balabdaoui
F.
Polakowski
M.
2005
Using Bayesian model averaging to calibrate forecast ensembles
.
Mon. Weather Rev.
133
(
5
),
1155
1174
.
http://dx.doi.org/10.1175/MWR2906.1
.
Sankarasubramanian
A.
Vogel
R. M.
2002
Annual hydroclimatology of the United States
.
Water Resour. Res.
38
(
6
),
1083
.
doi:10.1029/2001WR000619
.
Schneider
M. P. A.
Yaşar
Y.
2016
Is inequality deadly and for whom? a Bayesian model averaging analysis
.
The Social Science J.
53
,
357
370
.
Shepard
D.
1968
A two-dimensional interpolation function for irregularly-spaced data
.
Proceedings of the 23rd ACM National Conference
,
New York, pp. 517–524
.
Soares
T. N.
Diniz-Filho
J. A. F.
Nabout
J. C.
Telles
M. P. C.
Terribile
L. C.
Chaves
L. J.
2015
Patterns of genetic variability in central and peripheral populations of Dipteryx alata(FaBaceae) in the Brazilian Cerrado
.
Plant Syst. Evol.
301
(
5
),
1315
1324
.
Sullivan
D. O.
Unwin
D. J.
2003
Geographical Information Analysis
.
John Wiley & Sons, Inc.
,
N.J
.
Teegavarapu
R. S. V.
Elshorbagy
A.
2005
Fuzzyset-based error measure for hydrologic modelevaluation
.
J. Hydroinf.
7
,
199
208
.
Thiessen
A. H.
1911
Precipitation averages for large areas
.
Mon. Weather Rev.
39
(
7
),
1082
1084
.
Thomas
H. A.
1981
Improved Methods for National Water Assessment: Final Report
.
U.S. Geol. Surv. Water Resour. Contract WR15249270
, p.
44
.
Vieux
B. E.
2001
Distributed Hydrologic Modeling Using GIS, Water Science and Technology Library
.
Kluwer Academic Publishers
.
Wagner
P. D.
Fiener
P.
Wilken
F.
Kumar
S.
Schneider
K.
2012
Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions
.
J Hydrol.
464–465
,
388
400
.
Wang
F.
2006
Quantitative Methods and Applications in GIS
.
CRC Press
,
Boca Raton, Florida
,
USA
.
Wong
K. W.
Wong
P. M.
Gedeon
T. D.
Fung
C. C.
2003
Rainfall prediction model using soft computing technique
.
Soft. Comput.
7
,
434
438
.
Wu
D.
Cui
Y. L.
Xie
X. H.
Luo
Y. F.
2019
Improvement and testing of SWAT for multi-source irrigation systems with paddy rice
.
J Hydrol.
568
,
1031
1041
.
https://doi.org/10.1016/j.jhydrol.2018.11.057
.
Yang
H. B.
Yang
D. W.
Lei
Z. D.
Sun
F. B.
2008
New analytical derivation of the mean annual water-energy balance equation
.
Water Resour. Res.
44
(
3
).
http://dx.doi.org/10.1029/2007WR006135
.
Zhang
X.
Srinivasan
R.
2009
GIS-based spatial precipitation estimation: a comparison of geostatistical approaches
.
J. Am. Water Resour. Assoc.
45
(
4
),
894
906
.
Zhang
J.
Taflanidis
A. A.
2019
Bayesian model averaging for Kriging regression structure selection
.
Probabilist. Eng. Mech.
56
,
58
70
.
https://doi.org/10.1016/j.probengmech.2019.02.002
.
Zhao
J.
Wang
D.
Yang
H.
Sivapalan
M.
2016
Unifying catchment water balance models for different time scales through the maximum entropy production principle
.
Water Resour. Res.
52
.
https://doi.org/10.1002/2016WR018977
.
Zimmerman
D.
Pavlik
C.
Ruggles
A.
Armstrong
M. P.
1999
An experimental comparison of ordinary and universal kriging and inverse distance weighting
.
Math. Geol.
31
(
4
),
375
390
.