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

Spatial interpolation methods . | Group . | Abbreviation . |
---|---|---|

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 methods . | Group . | Abbreviation . |
---|---|---|

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

*p*at the unsampled spatial location is equal to the closest available rainfall,

_{u}*p*. If

_{i}*λ*is the distance between locations

_{u,i}*u*and

*i*, the weight of

*λ*is 1 for the location

_{u,i}*i*corresponding to the minimum distance,

*d*, and zero for all other locations (Thiessen 1911; Nikolopoulos

_{u,i}*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,

*λ*, are also related to the distance,

_{u,i}*d*, from the sampled point,

_{u,I}*p*, to the estimated point,

_{i}*p*

_{u}, as given by Equation (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).

*U*and

*V*are geographic coordinates,

*p*is a single observation of the mapped variable at coordinates (

_{i,j}*U*,

_{i}*V*). The

_{j}*p*at an unmeasured spatial location can be estimated by Equation (2) (Wong

_{u}*et al.*2003):where

*λ*is the polynomial coefficient and

_{i,j}*N*is the degree of the polynomial that represents the best surface of the precipitation pattern.

*λ*and

_{i,j}*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).

*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):where

*γ*(

*x*

_{i}*−*

*x*) represents the value of the semivariogram function for the distance between points

_{j}*x*and

_{i}*x*.

_{j}*γ*(

*x*

_{j}*−*

*x*) is the value of the distance between

_{u}*x*and the estimated location

_{j}*x*, and

_{u}*ϕ*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):

*C*

_{0},

*C*

_{0}+

*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:where the subscripts ‘1’ and ‘2’ refer to the primary and secondary networks, respectively, and

*N*

_{1}and

*N*

_{2}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).

*p*, the corresponding target, or training data,

*p*

_{T}, and an ensemble of all spatial interpolation method (i.e., TSA, NN, IDW, OK, and CK) simulations {

*f*

_{1},

*f*

_{2},

*f*

_{3},

*f*

_{4},

*f*

_{5}} of the precipitation. According to the law of total probability (Bulmer 1979), the probability density function of

*p*can be expressed as:where Pr(

*p*|

*f*) is the probability prediction given by the simulation of

_{k}*f*alone, and Pr(

_{k}*f*|

_{k}*p*) is the posterior probability of the method prediction, given the data set,

_{T}*p*. Pr(

_{T}*f*|

_{k}*p*) is also identified as a fractional statistical weight,

_{T}*w*, of the

_{k}*k*th spatial interpolation method, whose magnitude reflects how well

*f*matches

_{k}*p*, and . In order to determine the prior Pr(

_{T}*p*|

*f*) for the parameters and the prior probability of each method, it is computationally convenient to assume that Pr(

_{k}*p*|

*f*) is a Gaussian distribution defined by a mean,

_{k}*μ*, and variance,

_{k}*σ*

_{k}

^{2}(Raftery

*et al.*2005; Duan & Phillips 2010).where

*g*(•) is the Gaussian probability density function, and

*θ*= [

_{k}*μ*

_{k}, σ_{k}

^{2}] is the parameter vector.

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

*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.where

*O*and are the observation and estimation at the

_{i}*i*th gauge station, respectively.

#### Performance assessment of hydrological models

*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):

## 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 × 10^{6} km^{2} (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.

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

N*
. | Daily . | Monthly . | Annual . |
---|---|---|---|

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*
. | Daily . | Monthly . | Annual . |
---|---|---|---|

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.

Parameter . | Daily . | Monthly . | Annual . |
---|---|---|---|

N | 4 | 6 | 6 |

α | 1.6 | 1.8 | 1.6 |

Parameter . | Daily . | Monthly . | Annual . |
---|---|---|---|

N | 4 | 6 | 6 |

α | 1.6 | 1.8 | 1.6 |

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.

Semivariogram models . | Statistical criteria . | OK . | CK . | ||||
---|---|---|---|---|---|---|---|

Daily . | Monthly . | Annual . | Daily . | Monthly . | Annual . | ||

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 models . | Statistical criteria . | OK . | CK . | ||||
---|---|---|---|---|---|---|---|

Daily . | Monthly . | Annual . | Daily . | Monthly . | Annual . | ||

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 |

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.

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

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.

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

. | SWAT (Daily) . | abcd (Monthly) . | Budyko (Annual) . | |||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

. | NSE . | PBIAS(%) . | NSE . | PBIAS(%) . | NSE . | PBIAS(%) . | ||||||||||||

Spatial precipitation interpolation scheme . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . |

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 scheme . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . | I . | II . | III . |

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(a_{2}) showed more concordance between the simulated and observed values than in the calibration period (Figure 9(a_{1})). 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).

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