The best alternative for estimating reference crop evapotranspiration in different sub-regions of mainland China

Reference crop evapotranspiration (ETo) is a critically important parameter for climatological, hydrological and agricultural management. The FAO56 Penman-Monteith (PM) equation has been recommended as the standardized ETo (ETo,s) equation, but it has a high requirements of climatic data. There is a practical need for finding a best alternative method to estimate ETo in the regions where full climatic data are lacking. A comprehensive comparison for the spatiotemporal variations, relative errors, standard deviations and Nash-Sutcliffe efficacy coefficients of monthly or annual ETo,s and ETo,i (i = 1, 2, …, 10) values estimated by 10 selected methods (i.e., Irmak et al., Makkink, Priestley-Taylor, Hargreaves-Samani, Droogers-Allen, Berti et al., Doorenbos-Pruitt, Wright and Valiantzas, respectively) using data at 552 sites over 1961–2013 in mainland China. The method proposed by Berti et al. (2014) was selected as the best alternative of FAO56-PM because it was simple in computation process, only utilized temperature data, had generally good accuracy in describing spatiotemporal characteristics of ETo,s in different sub-regions and mainland China, and correlated linearly to the FAO56-PM method very well. The parameters of the linear correlations between ETo of the two methods are calibrated for each site with the smallest determination of coefficient being 0.87.

explicitly incorporate physiological and aerodynamic parameters 3,5,23 . The Penman-Monteith (PM) equation was selected as the standard ET o estimation method by the Food and Agriculture Organization (FAO) of the United Nations because it closely approximates ET o at the locations evaluated 5,24,25 . Afterwards the FAO56-PM equation was world widely applied for the validation of the other equations in absence of experimental measurements 26 .
In recent years, variations of FAO56-PM-based ET o (ET o,s ) has been extensively investigated since the FAO56-PM was recommended as a standard ET o estimation method [27][28][29][30][31] . The ET o,s variations were also analyzed in partial or entire mainland China (EMC) concerning different application objectives [32][33][34][35][36][37][38][39][40][41][42][43] . Meanwhile, the evaluation of the FAO56-PM method has also been widely conducted by comparing with different ET o estimation methods. The evaluation research mainly focused on answering which method could be an alternative for FAO56-PM either the input data were full, limited or missing 9,[44][45][46] . It is known that the FAO56-PM equation requires a large data input for estimating ET o , including the geological variables such as elevation and latitude, and the meteorological variables such as minimum air temperature (T min ), average air temperature (T ave ), maximum temperature (T max ), wind speed, relative humidity (RH) and sunshine hour (n). The high data demand of the FAO56-PM method realized its overall high accuracy, but restricted its application in some data-lacking regions. In the regions where the observed long-term meteorological data are difficult to obtain, the FAO56-PM method is not the best choice. To solve this problem, ET o estimation methods with a lower data requirement and a simpler computation process are preferentially applied.
Although ET p and ET o are not equivalent terms, both provide estimates of atmospheric evaporative demand. In the previous studies, there are different understandings about the relationship between ET p and ET o . Several researchers differ the two items strictly 32,47,48 . For instance, FAO56-PM ET o equation is considered a PM ET p equation for specific reference conditions 47 . For strictly utilization of the items, Allen et al. 5 strongly discourage the use of ET p for ET o estimation concerned about the ambiguities in their definitions. A few researchers consider ET o is a kind of ET p 49 or it is reference values of ET p for a uniform grass reference surface 50 . Noticeably, a lot of researchers look upon ET p and ET o as identical concepts and share similar equations for their estimations [51][52][53][54][55] . Usually, a climatologist or meteorologist and a hydrologist use the term "potential", whereas an irrigation scientist uses the term "reference crop", although the estimation equation could be same. Even some equation-proposers potentially identified the two items, such as Hargreaves and Samani 16, 17 adopted "potential" while Hargreaves and Samani 16 used the term "reference crop". Ambiguity between ET o and ET p was expected to be reduced by more extensive definition of ET p as potential crop evapotranspiration or by using one of the ET o definitions 56 . Take Thornthwaite 1 for another example, this method was originally proposed to estimate ET p , but was also applied for estimating ET o in different cases 57,58 . Therefore, although there are differences between ET p and ET o , there is close relationship between them and their estimations could be quantitatively linked.
China has a total land area of 9,597,000 km 2 and is the third largest country in the world. It has a complicated geomorphology which contains different water bodies, glaciers, frozen soils, deserts, basins, mountains, farmland, and forests. The elevations are general lower and lower from the west to the east, shaping a so called "3-level-catena" landform 59 . The weather stations in China distributed non-uniformly, there are more weather stations in the eastern China, but less in the western regions, especially less on Qinghai-Tibet Plateau. Neither is the distribution of the sites even, nor are the observed climatic elements same for different stations. The total sites available for air temperature and precipitation data are as large as 2474 60

Data and Methodology
Data. Geographical and weather data from 552 National Meteorological Observatory stations in EMC were collected from the China Meteorological Administration. The data contained both the daily and monthly timescales. The weather data included T max , T min , T ave , U 10 wind speed at 10 m, RH, and n. The data duration was 1961-2013. The elevations of the selected sites covered a large range in EMC (Fig. 1a). To obtain more accurate ET o estimation, the 48 sites reported by Chen et al. 62 were used as the radiation correction station. Meteorological station (marked with blue circle) and radiation calibration station (marked with red triangle) and they were set as the centers of the Thiessen polygons to find the other sites which would use same parameters with them for estimating radiation (Fig. 1b). The EMC is divided into seven sub-regions 63 considering the differences in topography and climate 64 (Fig. 1c) where G is soil heat flux (MJ m −2 month −1 ), T ave is mean air temperature at 2 m (°C), , U 2 and U 10 are wind speed at 2 and 10 m (m s −1 ), respectively, U 2 = 0.75 U 10 , e s is saturation vapor pressure (kpa), e a is actual vapor pressure (kpa), e s -e a is saturation vapor pressure deficit (kpa), Δ is slope of vapor pressure curve (kpa °C −1 ), γ is psychrometric constant (kpa °C −1 ), and R n is net radiation (MJ m −2 month −1 ). Monthly G is estimated by: where subscripts K + 1, K and K − 1 are order of month, respectively. Annual ET o,s is cumulated from the values of 12 months. R n is calculated by:  Performance evaluation of the 10 selected methods. Relative error (RE), standard deviation (θ) and Nash-Sutcliffe efficacy coefficient (NSE) 71 are used to assess the performances of monthly ET o,i : When NSE is close to 1, the quality of the method for estimating ET o,i is good with high reliability. When NSE is close to 0, ET o,i has an close mean value with ET o,s with an overall reliable estimation, but the errors of the estimation processes are large; when NSE is much less than 0, the estimation is not reliable.   6 Berti et al. 68 MHS_2  where r jj is sample autocorrelation coefficient of ET o,L at a lag jj. For denoting significance of a trend, when jj = 0, Z * equals to Z; while as jj > 0, the MMK statistic Z * is utilized.
Variation coefficient. The variability of series ET o,L is quantified with a coefficient of variation (C v ), calculated with the following equation (Nielsen and Bouma) 75 : where θ and ET L o, are standard deviation and multi-year mean ET o,L series, respectively. Variability levels are classified by C v ≤ 0.1, 0.1 < C v < 1.0 and C v ≥ 1.0 as weak, moderate or strong one, respectively.
Spatial distributions of the climatic variables, ET o,s , ET o,i and the other studied parameters are mapped by the Kriging interpolation method in ArcGIS 10.2 software.

Results
Spatial distribution of climatic variables. The spatial distribution of ET o are closely related to that of the related meteorological elements. Figure 2 illustrates the distribution of multi-year mean T ave , T min , T max , n, U 2 and RH. The distribution of T min , T max , T ave were generally similar, with high values in sub-regions V and VII but lower  Table 1. Detail description of the equations is in the section "Data and methodology". Figure 3 shows the spatial distribution of multi-year mean monthly ET o,s and trend test results of annual ET o,s series for each site. The site number for different trends is presented in Table 2. In Fig. 3a, the ET o,s values were higher in the sub-regions I, II, VI and VII than the other 3 sub-regions, ranging from 49 to 108 mm. ET o,s had a moderate spatial variability with a coefficient of variation (C v ) being 0.15. In general, ET o,s in western China (high elevations) were larger than in eastern and middle China (low elevations). In Fig. 3b and Table 2, more sites (339) had decrease trends in ET o,s than increase trends (213 sites), and the trends at more sites were insignificant. The sites which had decrease and increase trends occupied 61.4% and 38.6% of the total, respectively. This indicated an overall decrease of ET o,s in China. The common occurrence of insignificance trends was induced by the removing of autocorrelation structures when using the modified nonparametric Mann-Kendall test (MKK) method. It's reasonable for the trend analysis. The sites with significant decrease (Sig. Dec) trends in ET o,s were mainly located in eastern China, while the sites with significant increase (Sig. Inc) trends in ET o,s were mainly located in middle China (i.e., sub-region VI).     (Figs 5 and 6). In Fig. 5, the variation patterns of monthly ET o,i and ET o,s were general with single peak (valley) around July (January or December), which were also the months that the largest (smallest) differences between ET o,i and ET o,s occurred. The differences between ET o,i and ET o,s curves was the largest for the sub-region I (northwestern China), and was the smallest for the sub-region VI (the Qinghai-Tibetan Plateau). The   68 ), which belonged to the temperature-based type, needed only temperature data, and was simple in computation. In general, both ET o,i and ET o,s curves were regional-, seasonal-and method-specific.

The Spatiotemporal variation of ET
In  Because the estimated ET o,i values were regional-specific, the RE values for ET o,i also showed differences in spatial distributions (Fig. 7). Ranges of RE for ET o,i varied. The range of absolute RE values for ET o,9 was the largest, followed by ET o, 10 . RE for most of ET o,i covered both negative and positive values, but RE range of ET o,1 , ET o,9 and ET o,10 covered only negative values. ET o,7 had the smallest RE range, which reflected that the FAO24 method was more accurate for estimating monthly ET o in EMC. The radiation-based Mak method had smaller RE ranges when compared to the empirical, temperature-based methods and another radiation-based method PT, but it generally underestimated ET o in most of the months and sites and only had local adaptability in sub-region VI, therefore this method shouldn't be the best alternative for ET o,s in different sub-regions and EMC. Considering the simpler temperature-based ET o type, the MHS_2 method had lower RE than the other temperature-based methods. In general, the spatial distribution of RE for different ET o,i differed at different locations. It revealed the differences in adaptability extents of the applied methods.
Generally consistent with the spatial distribution, the temporal variations of RE for monthly ET o,i were also method and sub-region-specific (Fig. 8). The largest (smallest) RE curves generally occurred for   Each method overestimated or underestimated ET o in different months or the whole year when compared to FAO56-PM, but in the temperature type, the MHS_2 method was found to be the closest to FAO56-PM considering. Although the MHS_2 method underestimated the ET 0,s by 24% in January, and 15% and 25% in November and December, but had a very low RE (2%) for the year when compared to the FAO56-PM method.
Standard deviation. The spatial distribution of multiyear mean monthly standard deviation (θ) for ET o,i are illustrated in Fig. S2. All of the ten ET o estimation methods showed larger θ values in the northern China (sub-regions I, II and III), although with different ranges of θ. There was a method order of ranges for θ, i.e., IRA < Val_1 < M ak < PT < MHS_2 < HS < FAO24 < MHS_1 < Val_2 < KPM. In general, a larger ranges of ET o,i corresponded to a larger ranges of θ, the IRA method had a smallest range of θ because it had a smaller range of ET o. In fact, this method largely underestimated ET o,s . The KPM method had a largest range of θ, which indicated the variation  Table 4. θ in sub-region I, VI and EMC were generally smaller than the other sub-regions for each month. The MHS_1 had the largest θ for all of the sub-regions and EMC. For EMC, the θ curves ranked in the method order of Nash-Sutcliffe efficiency coefficient. The spatial (temporal) distribution of multiyear mean monthly Nash-Sutcliffe efficiency coefficients (NSEs) for ET o,i are illustrated in Figs S4 and S5. In Fig. S4, the ranges of NSE ranked in a method order of FAO24 < Mak < KPM < MHS_2 < PT < HS < MHS_1 < IRA < Val_2 < Val_1. The FAO24 and KPM, as analyzed above, were both combination based ET o methods, although both had smaller NSE ranges, their equations had higher demand of climatic variables. The Mak method had a smaller range of NSE than MHS_2 (between −9.32 and 0.35), it performed better than MHS_2 in sub-regions IV and VI, but it needed addition shortwave radiation (or sunshine hour) when estimating ET o . The NSE of MHS_2 method ranged between −16 and 0.20, it performed well for most sub-regions except VI and VII. From climatic variable demand aspect, the MHS_2 best met a simple equation standard than the other equations, with general good performance. In Fig. S5, the Val_2, Val_1, IRA, HS and Mak were excluded from the ten ET o methods because theire NSE values in each month and each sub-region were generally smaller than 0. Among the left 5 methods, similar to the RE values, the NSE of FAO24 and KPM methods were better, followed by the MHS_2 method, also indicating MHS_2's better performance for the temporal variations of monthly ET o in the methods with less climatic data demand.
NSE of the monthly and annual ET o,i using the selected 10 methods for EMC are presented in Table 5. Except ET o,7 which was estimated by the combination-based FAO24 method, there were 0, 0, 1, 2, 1, 4, 3, 0 and     68 in the non-combination type. For the whole year, only NSEs of the FAO24 and MHS_2 methods were larger than 0, which showed the feasibility of the two methods. But for a best alternative, the combination based FAO24 was not suitable. Therefore, through a comprehensive comparison of spatiotemporal variations from the ten selected methods by relative error, standard deviation and Nash-Sutcliffe efficiency coefficient, the MHS_2 method was preliminary selected as a better one for an alternative of ET o,s with its equation simplicity, least data demand and better performance.

Scatter plots of monthly ET o,i vs. ET o,s . Although the spatiotemporal distribution of multi-year mean
ET o,i were analyzed and the performances of all the methods were compared for each sub-region, direct comparison between monthly ET o,i and ET o,s are still necessary, in order that if the required full climatic data for estimating ET o,s are lacking, an relatively accurate alternative method could be selected out from the 10 candidate methods using less weather data. The scatter plots of monthly ET o,i with ET o,s for different sub-regions and EMC are illustrated in Fig. 9. In general, ET o,i deviated more with ET o,s in July, but less for December and January. In all of 7 sub-regions and EMC, Because temperature data are easier with less cost to observe, and ET o,6 estimated by the MHS_2 method was not only simpler, highly correlated with ET o,s in each month and most sub-regions, but also had generally good similarity in spatiotemporal distribution with ET o,s . Considering both good performance and the correlation with   where A and B are numerically equal to 1/a and −b/a, respectively. Equation 13 is also a calibration between ET o, 6 and ET o,s . For easier application of Eq. 13, values A and B for the 552 sites in China were validated. Figure 10 indicates the spatial distribution of A, B, and correlation coefficient (R). Values of A decreased from 1.32 to 0.67 from northwest to southwest and eastern China. B values were the largest in sub-region VI, followed by sub-regions II, III, IV, I, V, and VII, respectively. Values of R ranged between 0.87 and 0.99, were larger than 0.95 in most of China, especially in north China. The general high R values confirmed the applicability of the best alternative MHS_2 method in China after accurate calibration.

Discussion
Under the global climate change, decreasing trends in ET o have been observed in different parts of the world 32, 76, 77 , including China 78 and most parts of China, e.g., the Haihe River basin 79 , the Huang-Huai-Hai Plain 80 , the northwest China including Xinjiang Uywer Autonomos Region 81 , southeast China, the Yangtze river basin 64 , etc. The increasing trends were found at most sites of the Qinghai-Tibetan Plateau 34 . The trends were also bi-directional in China. This study revealed that annual ET o,s for 61.4% of the study sites had decreasing trends, of them, 9% of the trends were significant. Our research agreed with the former research in the general decreasing trends of ET o,s for EMC, but in the meantime, there were also differences between this research and the previous.. The differences were caused by the changes in the study period, the data source, the ET o,s estimation methods, the site number applied, and the research aims. For example, Wang et al. 51 also applied 4189-grid data during 1961-2013 in EMC to estimate ET o and identified the contribution of climatic variables to ET o variability. They revealed that annual ET o decreased with a mean rate of 6.84 mm/decade, and the sites with significant increase trends mainly distributed in the Qinghai-Tibetan Plateau. This research also reported general increasing trends in the same region, i.e. sub-region VI.
The most precise ET o estimation method varied for different regions. The frequently-used methods are the FAO56-PM, HS, and pan measurement etc., these methods have been applied to partial of China or EMC 58, 78, 82 . Xu et al. 82 applied 5 meteorological stations during 1999-2007 in arid-zone of China (i.e., sub-region I, VI of this research) and selected the HS method as the best alternative of ET o,s . This research selected the MHS_2 as a best alternative of ET o,s for different sub-regions and EMC, because it not only had a general high accuracy but also used only temperature data which were easy to observe or collected, even for the sites where the other climatic data were lacking. Moreover, this research also provided the spatial distributions of the calibrated parameters of  the MHS_2 method as the best alternative of ET o,s for different sub-regions and EMC, which were very useful for researchers to apply the calibrated MHS_2 method in China. The MHS_2 method overestimated ET o in the sub-regions V and VII in the high temperature section of EMC (Fig. 7f). RE reached 20% especially in March, April, May and June (Figs 8e,g and 9). Both sub-regions are humid and sub-humid climatic zones of EMC. This reflected the disadvantages of MHS_2 which only applied temperature data for estimating ET o . When temperature is high, ET o,6 obtained with the MHS_2 method could be high but ET o,s may not be as high as it considering also wind speed, relative humidity and sunshine hour. Under the overestimation conditions, the relationship between ET o, 6 and ET o,s should be re-calibrated for March, April, May and June. The re-calibrated parameters A, B and R 2 in March, April, May and June for the two sub-regions are presented in Table 7.  68 . The MHS_2 method utilized only temperature data, was simple in computation procedure when compared to the other 9 ET o estimation methods, and was highly correlated with ET o,s . It was a best alternative for ET o,s when climatic data were lacking. After accurate validation for the MHS_2 method using equation 13, the calibrated parameters of A and B for each site, sub-region and EMC were obtained. This research is an important contribution to ET o estimation method in China when the high requirements of climatic data could not be met.   Table 7. The re-calibrated parameters A, B and R 2 in March, April, May and June for the V and VII sub-regions using Equation 13.