#### Abstract

This study presented a method to estimate areal mean rainfall (AMR) using a Biased Sentinel Hospital Based Area Disease Estimation (B-SHADE) model, together with biased rain gauge observations and Tropical Rainfall Measuring Mission (TRMM) data, for remote areas with a sparse and uneven distribution of rain gauges. Based on the B-SHADE model, the best linear unbiased estimation of AMR could be obtained. A case study was conducted for the Three-River Headwaters region in the Tibetan Plateau of China, and its performance was compared with traditional methods. The results indicated that B-SHADE obtained the least estimation biases, with a mean error and root mean square error of −0.63 and 3.48 mm, respectively. For the traditional methods including arithmetic average, Thiessen polygon, and ordinary kriging, the mean errors were 7.11, −1.43, and 2.89 mm, which were up to 1027.1%, 127.0%, and 358.3%, respectively, greater than for the B-SHADE model. The root mean square errors were 10.31, 4.02, and 6.27 mm, which were up to 196.1%, 15.5%, and 80.0%, respectively, higher than for the B-SHADE model. The proposed technique can be used to extend the AMR record to the presatellite observation period, when only the gauge data are available.

#### 1. Introduction

Accurate estimation of areal mean rainfall (AMR) is necessary for water balance calculations, climate studies, flood modeling, and runoff forecasting [1]. Currently, rainfall data are derived mainly from three principal sources: (1) rain gauges; (2) weather radar; (3) and satellites [2–5]. However, the series of rainfall data based on weather radar and satellites are short, and therefore, to obtain the longest data series possible, AMR is mainly determined using rain gauge data.

Methods commonly used for the estimation of AMR can be classified into two categories. The first involves the calculation of the average or weighted average of rain gauge data, for example, using an arithmetic average method or Thiessen polygon technique [1], and the second requires the acquisition of the average precipitation for all pixels of a grid map using spatial interpolation methods such as one of the kriging techniques [5].

The arithmetic average method assumes that the rainfall field is homogeneous and that the rain gauge observations are independent; thus, it gives equal weight to all rain gauges [1, 6]. The Thiessen polygon technique sets a fixed weighting factor for the areal coverage of a station. It also assumes that the rainfall field is homogeneous within the area of coverage of each station [1, 7]. These assumptions are really only reasonable in primarily low-lying or flat terrain and/or where rainfall events are typically extensive and homogeneous (e.g., widespread frontal storms) [5].

Kriging methods (e.g., ordinary kriging, universal kriging, and cokriging) are known as geostatistical approaches or multivariate geostatistical techniques with spatial correlation patterns [6, 8–11]. Previous case studies have confirmed their effectiveness in the spatial interpolation of rainfall within mountainous terrain [11–14]. However, their estimation accuracy is influenced considerably by the density and distribution of the rain gauge network [15–17]. Generally, when the number of stations is sparse and/or their distribution is uneven, kriging methods may obtain biased estimations when lacking additional information from covariates (e.g., elevation).

Covariates, for example, latitude, longitude, and elevation are commonly used in kriging techniques [9–11]. They are also used in other interpolation methods, for example, regression techniques [18]. Many published high resolution precipitation climatologies, such as Parameter Regression on Independent Slopes Model (PRISM) and WorldClim, were developed by integrating explanation variables in their algorithms [18, 19]. However, uncertainty may also be the highest in mountainous and in poorly sampled areas [19]. For example, in a simple linear climate-elevation regression technique [18], available rain gauges often do not span the whole range of elevations in a region, especially in mountainous areas where most of the stations are located in river valleys or plains, which may introduce much uncertainty for spatial extrapolation. Satellite products, such as radar and Tropical Rainfall Measuring Mission (TRMM) data, are increasingly used as covariates in interpolation methods [20–23]. Satellite data are very promising, as they provide detailed spatial distribution information of rainfall, and can be transferred to different regions, especially in data scarce areas [23, 24].

In western China, difficulty of access and low population densities result in poor spatial coverage of rain gauges, particularly in the Gobi Desert or high mountains and plateaus. Thus, the distribution of rain gauge stations is sparse and uneven. Most gauges are located in river basins or valleys; however, rainfall varies considerably with elevation within mountainous areas [25]. Therefore, the estimated AMR will be biased when using the arithmetic average, Thiessen polygon, or kriging techniques. From the viewpoint of spatial sampling and statistical inference, the problem when estimating AMR using any of the abovementioned methods is that the biased samples cannot indicate the total population precisely. Thus, a proper estimator that is capable of obtaining an unbiased estimation of the population from the biased samples is required.

The Biased Sentinel Hospital based Area Disease Estimation (B-SHADE) model, which was developed to correct biased samples, held the ability to obtain the best linear unbiased estimation of the total (or mean) population [26–29]. It considers both the correlation between each sample and the population and the dependence between each pair of samples based on historical satellite data. Thus, the accurate AMR may be estimated using biased rain gauge observations in combination with satellite data based on this model. A case study was implemented by comparison with traditional methods including the arithmetic average, Thiessen polygon, and ordinary kriging. Such a technique can improve the accuracies of AMR estimation in remote areas where the rain gauge distribution is sparse and uneven and also extend the AMR record to the presatellite observation period, when only the gauge data are available.

#### 2. Study Area and Data

##### 2.1. Study Area

The Three-River Headwaters region was selected as the study area. This is located in the middle of the Tibetan Plateau, south of Qinghai Province in China, and it encompasses an area of about 360,000 km^{2}. It belongs to the plateau continental climate region. From northwest to southeast, the annual average temperature increases from −5.6 to 3.8°C; annual average precipitation increases from 262.2 to 772.8 mm, and the climatic zones change from semiarid to humid subtropical [30]. In addition, there are many huge mountains in the northwest of the region, where glaciers cover some areas [31].

There are only 14 available meteorological stations within this area, of which just two are located in the west. Most of the stations are located in river valleys (Figure 1), and their elevations range from 3323.2 to 4612.2 m a.s.l.; there are no stations higher than 5000 m a.s.l. (Table 1). Such a sparse and uneven distribution of meteorological stations is typical of western China (Figure 1).

##### 2.2. Data

###### 2.2.1. Rain Gauge Observations

Rain gauge data were downloaded from the China Meteorological Data Sharing Service System (http://cdc.nmic.cn). For the study period (January 1998 to December 2011), complete monthly rainfall datasets were acquired for 649 stations (Figure 1), 14 of which were located within the Three-River Headwaters region. The annual mean precipitation within the study area from 1998 to 2011 was between 332.6 and 728.2 mm (Table 1).

###### 2.2.2. TRMM Rainfall Estimates

For this study, we choose TRMM 3B43 (Version 7) monthly precipitation data from 1998 to 2011. This dataset was acquired from the National Aeronautics and Space Administration (NASA) website (http://trmm.gsfc.nasa.gov/). It has spatial resolution of [32, 33], which is selected as the research scale of this study; namely, the bias and comparison statistics used in the development, application, and validation of the techniques are calculated at the TRMM resolution of 0.25 deg lat/lon grid. It has been validated that the Version 7 TRMM 3B43 data have better accuracy than previous versions [34, 35]. Figure 2 shows the spatial distribution of TRMM-estimated annual mean precipitation from 1998 to 2011 in the Three-River Headwaters region. The annual mean precipitation for 1998–2011 ranged from 184.8 to 788.6 mm for TRMM data (Figure 2), while it ranged from 332.6 to 728.2 mm for rain gauges (Table 1); that is, precipitation range of TRMM is up to 52.6% higher than that from rain gauge observations.

###### 2.2.3. Relationship between TRMM Data and Rain Gauge Records

The TRMM data indicated high accuracy with respect to individual rain gauge observations. A scatter diagram of monthly TRMM-estimated rainfall against rain gauge observations from January 1998 to December 2011 demonstrates good linear correlation (Figure 3). Setting rain gauge precipitation as the independent variable and TRMM-estimated rainfall as the dependent variable, the linear regression through the origin was performed using SPSS 19 software. The results show that the slope of their linear regression function is 0.983 and the adjusted determination coefficient () is 0.946 (, ).

For each station, the TRMM-estimated rainfall showed a good linear relationship with the rain gauge observations. For each of the 14 stations within the study area, the linear regression through the origin was also performed between the monthly TRMM-estimated rainfall and rain gauge observations from 1998 and 2011. It is shown in Table 2 that the slopes of the regression functions range from 0.809 to 1.141; the adjusted determination coefficients () are all >0.92, and all of the regression functions passed the significance test at the 0.01 level.

#### 3. Methodology

##### 3.1. The Principle of the B-SHADE Model

The B-SHADE model is a model-assisted and data-driven method used for deviation correction of population estimation using biased samples [26]. Hu et al. [28] extended the original technique from total population estimation to mean population evaluation. They also clarified that historical unbiased samples (instead of the historical population) were effective when estimating the population mean. To achieve an unbiased estimation of a population using biased samples, the model takes into account both the “vertical” correlation and the “horizontal” dependence between the samples and the population. The “vertical” relevance between each sample and the population, which is represented by their ratio, is used to correct the sampling bias. The “horizontal” interaction between each pair of samples, which is the covariance between them, is used to increase the effectiveness of the sample size and reduce estimation variance [28]. Moreover, to guarantee a best linear unbiased estimation, each sample is allocated a unique weight. Then, the mathematical expectation of the population estimation should be equal to the true value of population, and the variance of the estimator should be minimized [26, 28].

As important input data for driving the B-SHADE model, a time series of historical unbiased sample set or population (historical data) is used for model construction and weights calculation. The estimated population value can be finally calculated by weighting sentinel samples. Being driven by the historical data, the B-SHADE model can be used to choose the best combination of sentinel sites according to the minimum estimation variance of the population, when the number of samples is given, or to calculate the smallest number of samples needed for a given estimation variance [28].

##### 3.2. Objective

The objective is to calculate the best linear unbiased estimation of AMR based on rain gauge observations. The discrete format true value of precipitation is expressed aswhere is the value of the th grid of TRMM gridded product and is the total number of grids within the area. The estimated AMR iswhere is the gauge (or alternatively corresponding TRMM grid) value of the th station, is the weight of the th station, and is the number of stations (). Accurate estimation of areal average precipitation in (2) should address two terms: the unbiased estimation conditionand the minimum estimation variancewhere is the statistical expectation and is the variance of . Considering (2), (3) can be written as

The relationship between the th station and AMR is expressed aswhere is the bias of the th station with respect to the areal mean . Considering (5), (6) can be written as

##### 3.3. Estimation of Weights

The objective can be finally converted into solving the combined equations (8) and (7) using the least-squares method:where is the covariance between the th station and the areal average value. is the number of unbiased historical samples of the population, and is the covariance between the th and th stations, both and are acquired based on historical TRMM-derived AMR data because of the high accuracy of the TRMM data, as discussed in Section 2.2.3, and is a Lagrange multiplier.

The weights can be calculated by

The theoretical variance can be calculated usingwhere is also evaluated using historical data [28].

After the weights of the rain gauges have been determined using (10), the AMR can be computed using (2). Even though the rain gauge observations are accurate, the discrepancy in scale between a gauge (point) and its corresponding TRMM grid (0.25 deg lat/lon) still exists, which will introduce potential uncertainties if the gauge values are used in calculating covariances. At the same time, as the TRMM data indicated high accuracy with respect to individual rain gauge observations, the TRMM values are used for covariance calculation to keep the consistency of data scale.

##### 3.4. Validation

To obtain the AMR estimation, both and in (8) should be determined based on historical TRMM-derived AMR data. In this study, we classified the monthly TRMM data into two intervals: January 1998 to December 2004 and January 2005 to December 2011. The data from the first interval were used to compute the weights that the B-SHADE model required, and the data from the second interval were used to compute the AMR as validation.

To evaluate the effectiveness of the AMR estimation using the B-SHADE model, three commonly used methods (arithmetic average, Thiessen polygon, and ordinary kriging) were chosen for comparison. The B-SHADE model was implemented using a software package which was developed by Hu et al. [28], and Thiessen polygon and ordinary kriging were carried out in ArcGIS 10.0 platform. The mean error (ME), root mean square error (RMSE), and mean absolute percent error (MAPE) were calculated using where and are the true and estimated AMR of the th month, respectively, and represents the number of months in the second 7-year study interval.

Two AMRs at two temporal scales were validated: the first was over a year and the second was over a season. Accuracies of all the four methods were also compared for optimum combinations of stations. For a given number of stations, the best combination of sentinel sites, which held the minimum estimation variance of the population (calculated based on historical TRMM data from 1998 to 2004), were selected from all possible combinations to compute the ME, RMSE, and MAPE for each technique.

#### 4. Results

##### 4.1. Multiyear Estimations of Monthly AMR

###### 4.1.1. Accuracies Based on Different Estimators

Table 3 shows that for different number of station combinations, the ME of the arithmetic average is always >3 mm, except when there are two stations. For the other methods, the MEs are smaller, especially for the B-SHADE model, for which the worst ME is only −2.30 mm. This indicates that the AMR is generally biased when estimated using the arithmetic average technique but that the other methods can correct the bias to different degrees. Among them, the B-SHADE model always has the least biased estimation.

Figure 4 shows the RMSEs of the arithmetic average, Thiessen polygon, and ordinary kriging fluctuate for the different numbers of stations. The RMSE of the arithmetic average is always the highest, followed by the kriging and Thiessen polygon techniques. The RMSE of the B-SHADE model is always the smallest in comparison with the other three methods. This indicates that the B-SHADE model obtains minimum estimation biases and the highest accuracies.

Table 3 and Figure 4 show that the accuracies of the AMRs estimated by the B-SHADE model are the highest of all four methods when all the rain gauge data are used. The ME and RMSE of the B-SHADE model are −0.63 and 3.48 mm, respectively. For arithmetic average, Thiessen polygon, and ordinary kriging, the mean errors were 7.11, −1.43, and 2.89 mm, which were up to 1027.1%, 127.0%, and 358.3%, respectively, greater than for the B-SHADE model; the root mean square errors were 10.31, 4.02, and 6.27 mm, which were up to 196.1%, 15.5%, and 80.0%, respectively, higher than for the B-SHADE model.

###### 4.1.2. Accuracy Change with an Increase of Rain Gauge Numbers

Accuracies of the estimations for different numbers of stations, by counting all the monthly precipitation from January 2005 to December 2011, show that when the number of stations increases, both the MEs and the RMSEs tend to decrease, except for the arithmetic average method (Table 3 and Figure 4). The ME and RMSE estimated by the arithmetic average model indicate no trend, varying along the averages of 5.4 and 9.2 mm, respectively. When the number of stations increases, the MEs and RMSEs of the arithmetic average in particular, but also of the Thiessen polygon and ordinary kriging methods, fluctuate severely (Table 3 and Figure 4). The RMSE of the arithmetic average method reaches its highest value of 14.39 mm when there are five stations and attains its lowest value of 5.71 mm when the number of stations is eight. Only the ME and RMSE of the B-SHADE method decrease gradually and reach their lowest values when all the rain gauge data are used. This indicates that the B-SHADE method has the most robust performance of the four methods.

##### 4.2. Estimation of Monthly AMR Variation with Season

###### 4.2.1. Accuracies Based on Different Estimators

Accuracies of estimations for different numbers of stations, by counting all the months every season from 2005 to 2011, show that both the ME and the RMSE estimated by the B-SHADE model generally indicate the lowest values, while those of the arithmetic average method have the largest values (Tables 4–7 and Figure 5). In summer, the accuracies of the estimations for the different numbers of stations using the B-SHADE model are the highest (Table 5 and Figure 5). In autumn, they are the highest except when only 2, 4, or 14 rain gauges are considered (Table 6 and Figure 5). In winter, the accuracies of the estimations using the B-SHADE model are similar to those based on the Thiessen polygon, but higher than the other two methods (Table 7 and Figure 5). In spring, the accuracies of the estimations using the B-SHADE model are similar to those based on ordinary kriging, lower than those using the Thiessen polygon method, and higher than those using the arithmetic average method (Table 4 and Figure 5).

**(a)**

**(b)**

**(c)**

**(d)**

###### 4.2.2. Accuracy Change with an Increase of Rain Gauge Numbers

Accuracies of the estimations for different numbers of stations, by counting all the months in one of the four seasons from 2005 to 2011, show that when the number of stations increases, both the MEs and RMSEs tend to decrease, except for the arithmetic average method (Figure 5). The ME and RMSE estimated by the arithmetic average model indicate no trend. The RMSEs of the arithmetic average vary along the averages of 5.57, 15.14, 8.68, and 0.84 mm in spring, summer, autumn, and winter, respectively, and the MEs vary along the averages of 6.0, 15, 9, and 0.8 mm in spring, summer, autumn, and winter, respectively.

When the number of stations increases, the MEs and RMSEs of the arithmetic average in particular, but also of the Thiessen polygon and ordinary kriging methods, fluctuate severely in every season. Only the MEs and RMSEs of the B-SHADE method decrease gradually and attain their lowest values when all the stations are considered. This is further evidence that the B-SHADE method has the most robust performance of the four methods at the season scale.

###### 4.2.3. Seasonal Variation of Accuracy Based on the B-SHADE Model

Figure 6 illustrates the MAPEs of the B-SHADE model in the different seasons. As the number of stations increases, the MAPE decreases gradually. When the number of stations is >10, the MAPE becomes relatively stable for all seasons. Among the four seasons, the MAPE displays the best accuracy in summer and poorest accuracy in winter. When the number of stations is >6, the accuracy in autumn is higher than in spring, whereas it is the opposite for smaller numbers of station combinations.

#### 5. Discussion

##### 5.1. Accuracy Improvement of AMR Estimation Combined with TRMM Data Based on B-SHADE Model

Based on the historical TRMM data, the population of the areal precipitation can be determined. Thus, through the historical TRMM data, the ratio of rainfall between each station and the population can be considered in calculating the weights for correcting the bias of the stations in the B-SHADE model. All the other methods take weighted averages to estimate the AMR; however, their weights are based simply on samples (rain gauges), which means that rainfall in areas where there is no rain gauge coverage cannot be reflected in the AMR.

Because the AMR bias correction using the B-SHADE model is based on both spatial and temporal distribution information of precipitation in the historical TRMM data, the B-SHADE model can be taken in AMR estimation in the pre-TRMM data period. If the spatial pattern information from the TRMM data can also be used by introducing covariates, the other three methods can perform better. However, the annual spatial patterns of precipitation in an area vary greatly. Thus, AMR estimation through introducing covariates cannot be extrapolated and be taken in the pre-TRMM data period.

##### 5.2. Robustness of AMR Estimation Based on B-SHADE Model

The B-SHADE model can select the best station combination that presents the smallest theoretical variance of precipitation from all the stations [26, 28]. The spatial pattern of rainfall based on the selected rain gauges can be a closer representation of the “real” one based on the TRMM data than other station combinations. Thus, the AMR will be the best linear unbiased estimation. Conversely, for the other three methods, there are no similar regimes for selecting the best station combination. In this study, the selected stations for the other three methods are the same as those based on the B-SHADE model. The bias of the AMR with the three methods cannot be corrected as it is in the B-SHADE model using historical TRMM data. If the bias of rainfall of the selected stations is high, the accuracy of the estimation will be low and vice versa. Thus, it is easy to understand that the accuracies of the estimated AMR can vary greatly with these three methods.

##### 5.3. Seasonal Variation of AMR Estimation Accuracy

Of all the four methods examined, the B-SHADE model does not illustrate the best accuracy in winter and spring, as it does in the other seasons. This might be due to the different spatial distribution of precipitation in these two seasons compared with that of the entire study period. In this study, the length of the TRMM time series is limited for each season (21 months), and it is insufficient to ensure stable estimations for the B-SHADE model. Therefore, we derived the weights for the AMR estimation on the season scale based on the monthly estimations of all four seasons over the first 7 years of the study period (84 months).

Strictly, there are only two seasons in the Three-River Headwaters region [36, 37]: the wet season from May to September and the dry season throughout the rest of the year. The wet season is relatively warm and brings about 80% of the annual precipitation. In contrast, the dry season is quite cold and brings little rainfall. Different atmospheric circulations and sources of moisture between the wet and dry seasons collectively lead to different distribution patterns of precipitation within the study area [38]. In the wet season, the southwest monsoon controls the weather within the study area, and moisture from the Indian Ocean is transported to the region by the southwesterly wind, which brings heavy rainfall. In contrast, during the dry season, the entire Tibetan Plateau is under the influence of a westerly circulation and, thus, long distances and topographic barriers restrict the transportation of moisture into the region, leading to little precipitation.

As the spatial distribution of precipitation in winter and spring (November to April) is different from that in summer and autumn (June to September) and the wet season brings about 80% of the annual precipitation, the weighting of the stations used in this study is largely decided by the spatial distribution of rainfall during this period, which conversely results in relatively poor performance of the B-SHADE model in winter and spring.

In addition, the poor performance of the B-SHADE model in winter and spring is also caused by the lower accuracies of TRMM satellite precipitation product in these seasons. During this period, snowfall is the primary precipitation form, and it is covered by snow and ice in the study area. However, as a representative satellite precipitation product, the TRMM data presents limited skills in capturing and quantifying precipitation in this period for its intrinsic inability of sensors [39–41]. The TRMM combines two different types of satellite sensors, that is, microwave (e.g., TRMM Microwave Imager) and infrared. Microwave sensors show weak ability to capture short rainfall events, because of the anomalous sampling of low earth-orbiting and the limited number of overpasses [39, 42]. Meanwhile, snow and ice surfaces may be interpreted as rainy clouds by microwave sensors, which may introduce errors in mountainous regions [41]. Besides, discrimination of rain/no rain category based on cloud top infrared temperature may also fail in mountainous regions [41, 42]. Therefore, when inaccurate TRMM data is used in the proposed model, the accuracy of the AMR estimation will be limited (see Section 5.4).

##### 5.4. Limitations of the AMR Estimation Based on the B-SHADE Model

As stated in Section 5.3, because the length of TRMM time series is limited for each season, we have to determine the weights for estimating the AMR based on the monthly estimations of all four seasons over the first 7 years of the study period. Furthermore, it has been validated that the B-SHADE model is effective in estimating reliable AMR on the multiyear scale. To estimate seasonal AMR for multiyear or annual scales, it is better to calculate the weights at the corresponding time scale, and historical data for a longer time series are necessary to obtain sufficient samples to ensure stable estimations by the B-SHADE model.

The B-SHADE model assumes that the rainfall features of an area are consistent within the time series. Because of the limitations of TRMM data, we used the data from January 1998 to December 2004 for model construction, and the data from January 2005 to December 2011 for validation. The validation period is very close to that of the model construction and may be not enough to validate the effectiveness of AMR estimation in the past decades. If the spatiotemporal patterns changed significantly during the study period, the estimated error might become larger. However, we believe the likelihood of substantial change in the spatiotemporal pattern within one century is not high. Therefore, as the records of most rainfall gauges are less than this time scale, historical AMR could be determined further using the B-SHADE method proposed in this study.

The performance of the B-SHADE model also relies on the accuracy of the satellite data, for it provides spatiotemporal distribution information of precipitation which is the basis of parameters calculation. Over- and underestimation of the TRMM data can both introduce errors to the AMR estimation. These errors which indicate inaccurate rainfall relationships between stations and the population can lead to incorrect ratios for bias correction and false weights for AMR estimation. Moreover, it might also mislead the determination of the best station combination. Therefore, prevalidation and potentially correction for a good linear relationship between TRMM data and rain gauge records are inevitable in the proposed method of this study.

#### 6. Conclusions

This study presented a method to estimate the AMR using the B-SHADE model with biased rain gauge observations and TRMM data, for remote areas where the distribution of rain gauges is sparse and uneven. Based on the B-SHADE model, the AMR estimation addresses two terms: the unbiased estimation condition between the rainfall from each gauge and the AMR, and the minimum estimation variances between two rain gauges. TRMM data were used to determine the population of the AMR for model construction and validation. The objectives were finally solved using the least-squares method.

A case study was conducted for the Three-River Headwaters region in the middle of the Tibetan Plateau in the south of Qinghai Province, China, which covers an area of about 360,000 km^{2}. The results indicated that the B-SHADE model obtained the least estimation biases, that is, the highest accuracy for the AMR estimation; when all the stations were used, the ME and RMSE were −0.63 and 3.48 mm, respectively. For arithmetic average, Thiessen polygon, and ordinary kriging, the mean errors were 7.11, −1.43, and 2.89 mm, which were up to 1027.1%, 127.0%, and 358.3%, respectively, greater than for the B-SHADE model; the root mean square errors were 10.31, 4.02, and 6.27 mm, which were up to 196.1%, 15.5%, and 80.0%, respectively, higher than for the B-SHADE model. When the number of stations increased, the MEs and RMSEs of the arithmetic average in particular, but also of the Thiessen polygon and ordinary kriging methods, fluctuated severely; only the B-SHADE method demonstrated robust performance.

Because the length of TRMM time series is limited for each season, the weights for the AMR estimation were determined based on the monthly estimations of all four seasons over the first 7 years of the study period. To estimate seasonal AMR over multiyear or annual scales, it would be better to calculate the weights at the corresponding time scale, and historical data for longer time series should be provided to obtain sufficient samples to ensure stable estimations by the B-SHADE model. The proposed technique can be to extend the AMR record to the presatellite observation period, when only the gauge data are available.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This study was supported by the National Program on Key Basic Research Project of China (nos. 2015CB954103 and 2015CB954101) and the National Key Technology Research and Development Program of the Ministry of Science and Technology of China (2012BAH33B01). The constructive criticisms and comments from the anonymous referees are also acknowledged.