Appropriateness of Potential Evapotranspiration Models for Climate Change Impact Analysis in Yarlung Zangbo River Basin, China

Evapotranspiration (ET) is an important element in the water and energy cycle. Potential evapotranspiration (PET) is an important measurement of ET. Its accuracy has significant influence on agricultural water management, irrigation planning, and hydrological modelling. However, whether current PET models are applicable under climate change or not, is still a question. In this study, five frequently used PET models were chosen, including one combination model (the FAO Penman-Monteith model, FAO-PM), two temperature-based models (the Blaney-Criddle and the Hargreaves models) and two radiation-based models (the Makkink and the Priestley-Taylor models), to estimate their appropriateness in the historical and future periods under climate change impact on the Yarlung Zangbo river basin, China. Bias correction methods were not only applied to the temperature output of Global Climate Models (GCMs), but also for radiation, humidity, and wind speed. It was demonstrated that the results from the Blaney-Criddle and Makkink models provided better agreement with the PET obtained by the FAO-PM model in the historical period. In the future period, monthly PET estimated by all five models show positive trends. The changes of PET under RCP8.5 are much higher than under RCP2.6. The radiation-based models show better appropriateness than the temperature-based models in the future, as the root mean square error (RMSE) value of the former models is almost half of the latter models. The radiation-based models are recommended for use to estimate PET under climate change in the Yarlung Zangbo river basin.


Introduction
With climate change and human activities, extreme hydrological events (such as drought and floods) happen more frequently and intensively, leading to a possible reduction of agriculture products and an increase of economic loss [1][2][3][4]. By 2050, human food demand will roughly double, due to population growth [5,6]. In response to this, improving the efficiency of use of irrigation water is useful [7,8]. Evapotranspiration (ET) plays a crucial role in both energy and water cycles, and its estimation accuracy has vital influence on irrigation planning, crop yield simulation, drought monitoring and forecasting [9][10][11][12][13][14].
The most usual expressions of ET are potential evapotranspiration (PET), actual evapotranspiration, and pan evaporation [15][16][17]. PET is the rate at which evapotranspiration occurs when the surface is well supplied with water [18]. Moreover, PET represents the comprehensive impact of various climate variables on the water use of vegetation [19,20]. Therefore, PET is often used in agricultural water

Study Area
The Yarlung Zangbo River, located on the north of the Himalaya Mountains, is the highest river over the world, with the elevation varying from 149 to 7100 m. The specific geographic position is between 29 • 9 N-31 • 9 N and 81 • 9 N-97 • 1 N (Figure 1). The basin area of the Yarlung Zangbo River Basin is about 246,000 km 2 and the length of this river is 2057 km (the fifth longest river in China). The dominant climate of this basin is a temperate semi-arid monsoon climate and the koppen-geiger climate classification of this territory is ET (E: polar, T: polar tundra) [55]. The moisture vapor of this basin comes mainly from the Indian Ocean [56]. The annual average precipitation ranges from 378 to 736 mm. The interannual variations of precipitation is quite tiny, but the distribution during the year is uneven. Precipitation throughout the year is mainly concentrated in July, August, and September, reaching 50-80% of annual precipitation. Moreover, the spatial distribution of precipitation differs greatly, with annual precipitation in the Gudi station from the upstream less than 300 mm and Baxika station from the downstream more than 4000 mm.
The annual mean air temperature is 5.3 to 7.6 • C, and it is likely to increase in the 21st century [57]. The studies of Wei and Fang (2013) [58] indicated that the temperature of the Yarlung Zangbo River Basin is rising. The reduction in the area of the glacier is the most significant manifestation of the increase in temperature [58]. Because of limited human activities [36], the Yarlung Zangbo River Basin is a desired study area to analyze the impacts of climate change on hydrological processes, including precipitation, runoff, and evapotranspiration.
Atmosphere 2019, 10, x FOR PEER REVIEW 3 of 25 climate change in the Yarlung Zangbo River Basin, and to investigate any uncertainty related to 18 global climate models under RCP2.6 and RCP8.5.

Study Area
The Yarlung Zangbo River, located on the north of the Himalaya Mountains, is the highest river over the world, with the elevation varying from 149 to 7100 m. The specific geographic position is between 29°9′ N-31°9′ N and 81°9′ N-97°1′ N ( Figure 1). The basin area of the Yarlung Zangbo River Basin is about 246,000 km 2 and the length of this river is 2057 km (the fifth longest river in China). The dominant climate of this basin is a temperate semi-arid monsoon climate and the koppen-geiger climate classification of this territory is ET (E: polar, T: polar tundra) [55]. The moisture vapor of this basin comes mainly from the Indian Ocean [56]. The annual average precipitation ranges from 378 to 736 mm. The interannual variations of precipitation is quite tiny, but the distribution during the year is uneven. Precipitation throughout the year is mainly concentrated in July, August, and September, reaching 50-80% of annual precipitation. Moreover, the spatial distribution of precipitation differs greatly, with annual precipitation in the Gudi station from the upstream less than 300 mm and Baxika station from the downstream more than 4000 mm.
The annual mean air temperature is 5.3 to 7.6 °C, and it is likely to increase in the 21st century [57]. The studies of Wei and Fang (2013) [58] indicated that the temperature of the Yarlung Zangbo River Basin is rising. The reduction in the area of the glacier is the most significant manifestation of the increase in temperature [58]. Because of limited human activities [36], the Yarlung Zangbo River Basin is a desired study area to analyze the impacts of climate change on hydrological processes, including precipitation, runoff, and evapotranspiration.

Data
The observation data used in this study were collected from 12 meteorological stations operated by the Tibet Meteorological Administration. The time period of data is 1971-2000, also used as the baseline period. A widely used meteorological baseline period is a 30-year "normal" period, as defined by the World Meteorological Organization (WMO) [59]. The main climatic variables include maximum and minimum air temperature, precipitation, wind speed, relative humidity, and sunshine hour at daily time step. Figure 1 shows the location of the 12 meteorological stations, and their elevation ranges from 2992 to 4900 m.
The simulated data by 18 GCMs on monthly time step were downloaded from the website http://cmip-pcmdi.llnl.gov/cmip5/. The detailed information of 18 GCMs are shown in Table 1. The

Data
The observation data used in this study were collected from 12 meteorological stations operated by the Tibet Meteorological Administration. The time period of data is 1971-2000, also used as the baseline period. A widely used meteorological baseline period is a 30-year "normal" period, as defined by the World Meteorological Organization (WMO) [59]. The main climatic variables include maximum and minimum air temperature, precipitation, wind speed, relative humidity, and sunshine hour at daily time step. Figure 1 shows the location of the 12 meteorological stations, and their elevation ranges from 2992 to 4900 m. The simulated data by 18 GCMs on monthly time step were downloaded from the website http://cmip-pcmdi.llnl.gov/cmip5/. The detailed information of 18 GCMs are shown in Table 1. The emission scenarios used in this study are RCP2.6 (low emission) and RCP8.5 (high emission). The future time period was chosen as 2041-2070.

Methodology Framework
The methodology framework of this study is shown in Figure 2. The first step is to locally calibrate parameters from four PET models for each station. The Penman-Monteith model is used here as a benchmark. The calibrated models will be used in the baseline and future periods. In the second step, bias correction is implemented for GCM output. The QM-Gaussian distribution method is used for maximum and minimum temperature, and the monthly correction factor (MCF) method is used for relative humidity, wind speed, and solar radiation. In the third step, the performance of the different PET models is compared in the baseline period, and changes of annual, seasonal, and monthly PET based on the different models are analyzed in the future period. Finally, the appropriateness of the PET models under climate change and uncertainty of global climate models are investigated.

PET Models
In this study, five widely used PET models were chosen, including one combination model (the FAO Penman-Monteith model), two temperature-based models (the Blaney-Criddle and the Hargreaves models) and two radiation-based models (the Makkink and the Priestley-Taylor models).

FAO Penman-Monteith Model
The Penman family of models is generally thought to be able to estimate PET accurately in most climates [11,29]. The FAO version of the Penman-Monteith model is recommended to calculate PET, owing to its accuracy and reality [27]. The high requirement of meteorological input is the main limitation of this model, thereby its use is often limited in data-sparse areas. The formula for the FAO Penman-Monteith model is:

PET Models
In this study, five widely used PET models were chosen, including one combination model (the FAO Penman-Monteith model), two temperature-based models (the Blaney-Criddle and the Hargreaves models) and two radiation-based models (the Makkink and the Priestley-Taylor models).

FAO Penman-Monteith Model
The Penman family of models is generally thought to be able to estimate PET accurately in most climates [11,29]. The FAO version of the Penman-Monteith model is recommended to calculate PET, owing to its accuracy and reality [27]. The high requirement of meteorological input is the main limitation of this model, thereby its use is often limited in data-sparse areas. The formula for the FAO Penman-Monteith model is: where PET is the potential evapotranspiration in mm/day; ∆ is the slope of the vapor pressure-temperature curve in kPa/ • C; R n is the net radiation at the crop surface in MJm −2 d −1 ; G is soil heat flux density at the soil surface in MJm −2 d −1 ; γ is the psychrometric constant in kPa/ • C; T is the mean daily air temperature at 2 m height in • C; u 2 is the wind speed at 2 m height in m/s; e s is the saturation vapor pressure at 2 m height in kPa and e a is the actual vapor pressure at 2 m height in kPa.

Blaney-Criddle Model
The Blaney-Criddle (1950) [24] model for estimating PET was originally proposed in the western United States. This model is also widely used in other regions of the world. The usual form of the Blaney-Criddle model is written as: where PET is the potential evapotranspiration in mm/day; T a is the mean air temperature in • C; p is the ratio of total daytime hours out of total daytime hours of the year; and k is the monthly consumptive use coefficient, relating to season, location, and vegetation type. The range of k is 0.5-1.2 for the growing period and an average value of 0.85 will be used as the initial value. The final value is calibrated locally.

Hargreaves Model
The Hargreaves model is a temperature-based model for estimating PET, and only temperature data and latitude are needed [60,61]. Owing to its low requirements of meteorological input, the Hargreaves model is widely used in many regions. The formula for the Hargreaves model is described as where PET is the potential evapotranspiration in mm/day; R A is the extraterrestrial radiation in mm/day; TD is the difference between maximum and minimum temperatures in • C; T a is the mean air temperature in • C; and k is an empirical coefficient. The initial value of k is 0.0023, and the final value is calibrated locally.

Makkink Model
Makkink (1957) [25] proposed a model to calculate grassland potential evapotranspiration under the cold climate of Holland. According to the experimental data of Holland, Hansen (1984) [62] modified the Makkink model. Recently, the Makkink model is widely used in western Europe, the USA, and other regions [63][64][65]. The formula for the Makkink model is as follows: where PET is the potential evapotranspiration in mm/day; ∆ is the slope of the vapor pressure-temperature curve in kPa/ • C; γ is the psychrometric constant in kPa/ • C; R s is the solar radiation in cal/(cm 2 ·d); λ is the latent heat in cal/g; and k is the empirical coefficient. The initial value of k is 0.7 and the final value can be calibrated locally.

Priestley-Taylor Model
Priestley and Taylor (1972) [66] proposed a simplified version of the Penman model for estimating PET from an extensive wet surface. In this model, the aerodynamic part was removed and the energy part was multiplied by a coefficient, α. The Priestley-Taylor model is as follows where PET is the potential evapotranspiration in mm/day; ∆ is the slope of the vapor pressure-temperature curve in kPa/ • C; γ is the psychrometric constant in kPa/ • C; R n is the net radiation at the crop surface in cal/(cm 2 ·d); λ is the latent heat in cal/g; and α is the empirical coefficient, depending on the wetting condition. The initial value of α is 1.26 and the final value can be calibrated locally.

Projection of Future PET and Performance Evaluation
Future potential evapotranspiration is obtained by the delta change method [42,67,68], i.e., computing the differences between baseline and future GCM-derived PET, and adding these changes to PET based on observation data in the baseline period. This method is able to directly avoid the possible effects of bias in climate variables from GCMs on PET simulations. The formulas are as follows: where PET GCM− f uture is the potential evapotranspiration calculated based on GCM simulations in the future period and PET GCM−baseline is the potential evapotranspiration calculated based on GCM simulations in the baseline period. The performance of the PET models was evaluated using the percentage of absolute bias (PBIAS) and root mean square error (RMSE), which are respectively defined as follows: where E i and R i refer to the reference (calculated by the FAO-PM model) and estimated values (calculated by other models) of PET at the ith time step, respectively; n refers to the number of time steps; and − R refers to the mean value of reference PET values. The values of perfect fit for performance evaluation are both 0 for PBIAS and RMSE.

Bias Correction Method
The output from GCMs or RCMs are generally biased, making correction procedures necessary before they are applied to regional climate impact analysis [44,45]. Bias correction methods identify possible bias between observed and simulated data and then adjust the simulated data. In this study, the quantile mapping (QM) method is used to do bias correction for maximum and minimum temperatures and the monthly correction factor (MCF) method is used for relative humidity, wind speed, and solar radiation.
The QM method is a non-parametric bias correction method proposed by Themessl et al., (2011) [53]. This method was successfully conducted in the bias correction for GCMs or RCMs-simulated precipitation in many studies [44,45]. Wilcke et al. (2013) [69] adopted the QM method to correct the bias of six climate variables from GCMs and results indicated that annual and monthly biases can largely be reduced for all variables. For precipitation, Gamma or double Gamma distribution is often used in the QM method. However, for temperature, a Gaussian distribution is usually assumed to match temperature better [38,67]. Hence, QM-Gaussian distribution is implemented for the bias correction of maximum and minimum temperature in this study.
In most cases, only simulated precipitation and temperature data are bias corrected. However, biases also occur in other variables, such as relative humidity, wind speed, and solar radiation [54]. These variables play vital roles on PET estimation [51,52]. The MCF method proposed by Haddeland  [54] is used here to correct the bias in other variables. The MCF method is performed at a monthly time step as follows: where V bc is the variable value after bias correction; V sim is the raw GCM output variable value; V obs (m) is the long-term mean monthly observed value of variable; and V sim (m) is the long-term mean monthly simulated value of variable. The long-term mean differences ( ) in the period 1971-2000 were used for baseline and future periods. The performance evaluation for the above bias correction methods is root mean square error (RMSE) [70]. The equation is listed as Equation (9).

Comparison of PET in the Baseline Period
The parameters from the four models were calibrated locally for each station. The Penman-Monteith model was used here as a benchmark. Root mean square error (RMSE) was used as the objective for parameter calibration of the four models shown in Section 3.2. The parameter corresponding to the smallest RMSE was regarded as the optimal parameter. Since the Blaney-Criddle model was originally developed for seasonal use [24], its parameter was then calibrated at seasonal time step. The other three models were calibrated at yearly time step. Table 2 shows the calibrated parameters in the four models for each station. As shown in Table 2, big differences exist between the values of the original parameter and the calibrated parameter.   Figure 3 shows a comparison of performance evaluation indices from annual PET before and after calibration. From this figure, it can be shown that the PBIAS and RMSE decrease largely after calibration, particularly for temperature-based models which display large biases before calibration. For instance, the performance evaluation of the Blaney-Criddle model at Longzi Station (No. 8) improved greatly after calibration, decreasing from 13.53% to 0.79% and from 0.42 to 0.096 mm/day for PBIAS and RMSE, respectively. This fact demonstrates that parameter calibration is indeed necessary. It was found out that the performance from the twelve stations was different, especially at Pulan Station (No. 10) where the four models revealed bad performance with the FAO-PM model. The reason was that wind speed is very important in estimating PET [51] but was not included in the four models, and the average value of wind speed at the Pulan station is much bigger than other stations. Generally, the performance evaluation indices (PBIAS and RMSE) indicate that annual PET estimated by the Blaney-Criddle and Makkink models are in better agreement with those from the FAO-PM model than Hargreaves and Priestley-Taylor models.
improved greatly after calibration, decreasing from 13.53% to 0.79% and from 0.42 to 0.096 mm/day for PBIAS and RMSE, respectively. This fact demonstrates that parameter calibration is indeed necessary. It was found out that the performance from the twelve stations was different, especially at Pulan Station (No. 10) where the four models revealed bad performance with the FAO-PM model. The reason was that wind speed is very important in estimating PET [51] but was not included in the four models, and the average value of wind speed at the Pulan station is much bigger than other stations. Generally, the performance evaluation indices (PBIAS and RMSE) indicate that annual PET estimated by the Blaney-Criddle and Makkink models are in better agreement with those from the FAO-PM model than Hargreaves and Priestley-Taylor models.        Taking the sensitivity of PET to climate variables into consideration, the reason for the higher deviation is probably because the four models do not consider wind speed, which has a significant impact on PET [51]. Moreover, the average value of wind speed at the Pulan station is much bigger than other stations, which leads to the higher deviation. In general, the Blaney-Criddle and Makkink models revealed better appropriateness than the Hargreaves and Priestley-Taylor models at a monthly time step. This conclusion is the same for the seasonal performance of the four models.  As presented in Figure 6, it can be found that spatial heterogeneity exists in the appropriateness of the four models for estimating monthly PET. In January, the Hargreaves and Priestley-Taylor models obtained lower PET values in the entire river basin, which is consistent with their performance in winter (DJF). In April, the Hargreaves model obtained lower PET values in the upstream and midstream of the river basin. In July, the Blaney-Criddle model showed the closest agreement with the FAO-PM model, while the Hargreaves model overestimated the PET values in the upstream and midstream of the river basin. In October, the spatial patterns based on the five PET models were rather similar, except the Makkink model. As presented in Figure 6, it can be found that spatial heterogeneity exists in the appropriateness of the four models for estimating monthly PET. In January, the Hargreaves and Priestley-Taylor models obtained lower PET values in the entire river basin, which is consistent with their performance in winter (DJF). In April, the Hargreaves model obtained lower PET values in the upstream and midstream of the river basin. In July, the Blaney-Criddle model showed the closest agreement with the FAO-PM model, while the Hargreaves model overestimated the PET values in the upstream and midstream of the river basin. In October, the spatial patterns based on the five PET models were rather similar, except the Makkink model.  Figure 7 shows the RMSE value between observation data with raw and corrected GCM maximum temperature and solar radiation in the baseline period. As presented in Figure 7, the RMSE values of two climate variables decrease significantly after bias correction for all GCMs. The range of RMSE is 2.3-17.7 °C for raw GCM maximum temperature data, while the range is 1.5-8.8 °C after bias correction. The ecdf of observed maximum temperature and corrected GCMs data also show close agreements for all GCM models (see Figure S1). This confirms the capability of the QM-Gaussian method for the bias correction of temperature data. For other climate variables (solar radiation, relative humidity, and wind speed), the QM method was actually adopted to correct bias of them, but the effect was unsatisfactory. Thus, the MCF method proposed by Haddeland et al., (2012) [54] was chosen for these climate variables. After the MCF method used, the values of RMSE from all GCM models were greatly reduced. The value of RMSE was the range of RMSE of 1.98-7.67 MJ/m 2 day for raw GCM solar radiation data, while this value was close to zero at most stations after bias correction. The bias correction of wind speed and relative humidity by the MCF method also yielded superior results. It is found that ecdf of corrected GCMs data are much closer to observation data for most GCM models (see Figure S2).  Figure 7 shows the RMSE value between observation data with raw and corrected GCM maximum temperature and solar radiation in the baseline period. As presented in Figure 7, the RMSE values of two climate variables decrease significantly after bias correction for all GCMs. The range of RMSE is 2.3-17.7 • C for raw GCM maximum temperature data, while the range is 1.5-8.8 • C after bias correction. The ecdf of observed maximum temperature and corrected GCMs data also show close agreements for all GCM models (see Figure S1). This confirms the capability of the QM-Gaussian method for the bias correction of temperature data. For other climate variables (solar radiation, relative humidity, and wind speed), the QM method was actually adopted to correct bias of them, but the effect was unsatisfactory. Thus, the MCF method proposed by Haddeland et al., (2012) [54] was chosen for these climate variables. After the MCF method used, the values of RMSE from all GCM models were greatly reduced. The value of RMSE was the range of RMSE of 1.98-7.67 MJ/m 2 day for raw GCM solar radiation data, while this value was close to zero at most stations after bias correction. The bias correction of wind speed and relative humidity by the MCF method also yielded superior results. It is found that ecdf of corrected GCMs data are much closer to observation data for most GCM models (see Figure S2).  Figure 8 presents the relative changes of future annual PET estimated by the five models at each station under two RCPs. It is easy to found that the relative changes of annual PET under RCP8.5 are much higher than RCP2.6. The Blaney-Criddle and Hargreaves models both project positive changes in PET at all GCMs. However, the Makkink and Priestley-Taylor models project negative future PET changes at a few GCMs, such as MIROC-ESM (No.14) and MIROC-ESM-CHEM (No.15). As shown in Figure 8, it can be observed that the difference among stations is much smaller than differences  Figure 8 presents the relative changes of future annual PET estimated by the five models at each station under two RCPs. It is easy to found that the relative changes of annual PET under RCP8.5 are much higher than RCP2.6. The Blaney-Criddle and Hargreaves models both project positive changes in PET at all GCMs. However, the Makkink and Priestley-Taylor models project negative future PET changes at a few GCMs, such as MIROC-ESM (No.14) and MIROC-ESM-CHEM (No.15). As shown in Figure 8, it can be observed that the difference among stations is much smaller than differences among GCMs. Generally, the change signs of future annual PET under two RCPs using the five models are consistent at most GCMs. In conclusion, the relative change of future annual PET under RCP2.6 estimated by the Priestley-Taylor model shows the closest agreement with those by the FAO-PM model, followed by Makkink, Hargreaves, and Blaney-Criddle models. Additionally, the performance of the five models under RCP8.5 is difficult to judge.   Figure 9 presents RMSE between seasonal PET change estimated by the FAO-PM model and the other four models in the future period under two RCPs. Similar to seasonal PET, RMSE under RCP8.5 is much higher than under RCP2.6. The appropriateness of each model performed differently with four seasons. The four models performed better in autumn (SON) and winter (DJF) than spring (MMA) and summer (JJA), which was contrary to that of baseline period. As presented in Figures 9  and 10 the sum RMSE values of four seasons from the radiation-based models is lower than those of the temperature-based models. Similar to the annual time step, the appropriateness of radiation-based models at seasonal time step was much better than that of temperature-based models in the future period. The other performance evaluation index PBIAS also confirmed this conclusion.

Changes in Annual, Seasonal and Monthly PET
Atmosphere 2019, 10, x FOR PEER REVIEW 15 of 25 Figure 9 presents RMSE between seasonal PET change estimated by the FAO-PM model and the other four models in the future period under two RCPs. Similar to seasonal PET, RMSE under RCP8.5 is much higher than under RCP2.6. The appropriateness of each model performed differently with four seasons. The four models performed better in autumn (SON) and winter (DJF) than spring (MMA) and summer (JJA), which was contrary to that of baseline period. As presented in Figure 9 and Figure 10 the sum RMSE values of four seasons from the radiation-based models is lower than those of the temperature-based models. Similar to the annual time step, the appropriateness of radiation-based models at seasonal time step was much better than that of temperature-based models in the future period. The other performance evaluation index PBIAS also confirmed this conclusion.  Figure 10 shows the relative changes of the average monthly PET values of 18 GCMs estimated by the five models. Obviously, the relative changes of PET calculated by the five models under RCP8.5 are much higher than RCP2.6, especially the Blaney-Criddle and Hargreaves models. This may be because these two models are temperature-based models and are, therefore, more sensitive to change of temperature in the future. As presented in Figure 10, the change signs of PET under two RCPs using the five models are consistent at the majority of months. The relative changes projected by temperature-based models from October to April are much larger than that of the FAO-PM model.  Figure 10 shows the relative changes of the average monthly PET values of 18 GCMs estimated by the five models. Obviously, the relative changes of PET calculated by the five models under RCP8.5 are much higher than RCP2.6, especially the Blaney-Criddle and Hargreaves models. This may be because these two models are temperature-based models and are, therefore, more sensitive to change of temperature in the future. As presented in Figure 10, the change signs of PET under two RCPs using the five models are consistent at the majority of months. The relative changes projected by temperature-based models from October to April are much larger than that of the FAO-PM model. The relative changes of monthly PET under two RCPs estimated by the Priestley-Taylor and Makkink models show closer agreement with those of the FAO-PM model than the Hargreaves and Blaney-Criddle models at most stations, especially at the Jiali and Gaize stations. Similar to annual and seasonal time steps, a conclusion can be drawn that the radiation-based models show better appropriateness with FAO-PM than the temperature-based models.

Uncertainty of Global Climate Models
To investigate the uncertainty associated with the 18 GCMs, boxplots of monthly areal PET relative change values estimated by the five models in the future period under two RCPs are presented in Figure 12. In this figure, the line in the middle of the box refers to the median value. The whiskers indicate the maximum and minimum values. The crosses outside the boxes represent outliers. As shown in Figure 12, it can be observed that the uncertainty under RCP8.5 was slightly higher than those under RCP2.6. For PET estimated by the FAO-PM model, higher uncertainty occurred in summer and winter. For the Blaney-Criddle model, relative changes of PET from June to October showed lower uncertainty than in other months, which was the same for the Hargreaves model. PET estimated by Makkink and Priestley-Taylor models illustrated similar characteristics of uncertainty, i.e., higher uncertainty occurs in summer and winter. This means that the uncertainty of PET from the same category reveals similar characteristics. The uncertainty from the temperaturebased models (Blaney-Criddle and Hargreaves) were the smallest, followed by the radiation-based models (Makkink and Priestley-Taylor) and the combined model (FAO-PM).

Uncertainty of Global Climate Models
To investigate the uncertainty associated with the 18 GCMs, boxplots of monthly areal PET relative change values estimated by the five models in the future period under two RCPs are presented in Figure 12. In this figure, the line in the middle of the box refers to the median value. The whiskers indicate the maximum and minimum values. The crosses outside the boxes represent outliers. As shown in Figure 12, it can be observed that the uncertainty under RCP8.5 was slightly higher than those under RCP2.6. For PET estimated by the FAO-PM model, higher uncertainty occurred in summer and winter. For the Blaney-Criddle model, relative changes of PET from June to October showed lower uncertainty than in other months, which was the same for the Hargreaves model. PET estimated by Makkink and Priestley-Taylor models illustrated similar characteristics of uncertainty, i.e., higher uncertainty occurs in summer and winter. This means that the uncertainty of PET from the same category reveals similar characteristics. The uncertainty from the temperature-based models (Blaney-Criddle and Hargreaves) were the smallest, followed by the radiation-based models (Makkink and Priestley-Taylor)

Discussion
To detect the appropriateness of PET models under climate change, this study evaluated the performance of five PET models in baseline and future periods, including one combination model (FAO-PM), two temperature-based models (Blaney-Criddle and Hargreaves), and two radiation-based models (Makkink and Priestley-Taylor). FAO-PM is used here as a benchmark to locally calibrate the parameters of the other four models. Local calibration for reference crop evapotranspiration models is also suggested from another study [19]. The results revealed that the Blaney-Criddle and Makkink models show better appropriateness in the baseline period. However, the appropriateness of PET models alters under climate change, i.e., the radiation-based models (Makkink and Priestley-Taylor) showed closer agreement with FAO-PM than the temperature-based models (Blaney-Criddle and Hargreaves). This confirms the conclusion from a previous study by Liu et al., (2017) [19] that radiation-based models perform better than temperature-based models for reference crop evapotranspiration at a semiarid site in China.
As shown in Figure 10, it can be concluded that future monthly PET estimated by temperature-based models is much higher than that obtained by FAO-PM, while radiation-based models show lower values. To illustration this phenomenon, we take Zedang Station as an example. Figure 13 presents the future average monthly change of the 18 GCMs for five climate variables at Zedang Station under two RCPs. Figure 14 shows the contribution percentage of climate variables to PET (estimated by FAO-PM) at Zedang Station based on the sensitivity analysis from another study by the authors [51]. As shown in Figure 14, PET is most sensitive to wind speed, followed by relative humidity, maximum temperature, solar radiation, and minimum temperature. Moreover, the contribution percentages vary at different months. Therefore, it can be figured out that the higher PET values calculated by the temperature-based models are mainly owing to the fact that these models only consider temperature and the significant increase of temperature in the future. The total contribution percentage of maximum and minimum temperatures ranges from 16% to 48% and a lower value occurs from May to August, demonstrating that the variations of temperature has limited impact on future PET estimation. Hence, the future PET values calculated by temperature-based models is higher than that obtained by the FAO-PM model, particularly from May to August. As presented in Figure 13, wind speed, relative humidity, and solar radiation have small decreases in the future period. The decrease of solar radiation and wind speed make PET decrease, and the decrease of relative humidity makes PET increase. Considering that the contribution percentage of solar radiation is small and its slight decrease, the change impact of solar radiation can be ignored in the future period. Nevertheless, mean temperature data is needed for the computation of variable ∆ (slope of the vapor pressure-temperature curve), which is used in the radiation-based and combination models. Thus, future PET estimated by the radiation-based models revealed a positive change trend, but the values were lower than that obtained by FAO-PM. For the other two variables, the decreased extent of wind speed and its contribution percentage were bigger than relative humidity, therefore, the impact of these two variables make future PET increase slightly. This interprets why PET calculated by FAO-PM is higher than that obtained by the radiation-based models. . Figure 14. Contribution percentage of climate variables to PET by FAO-PM at Zedang Station.
In this study, the QM-Gaussian method was used to correct the bias in maximum and minimum temperatures from GCMs and the results are very satisfactory. Although bias correction for temperature and precipitation is often made in other studies, bias correction for other climate variables is rarely conducted. Several previous studies have developed new tools or used existing approaches to bias correct for other climate variables, such as the distribution-based scaling (DBS) approach for relative humidity and wind speed [71], the QM method for all climate variables [44], and the MCF method for solar radiation, relative humidity, and wind speed [54]. In this study, the  In this study, the QM-Gaussian method was used to correct the bias in maximum and minimum temperatures from GCMs and the results are very satisfactory. Although bias correction for temperature and precipitation is often made in other studies, bias correction for other climate variables is rarely conducted. Several previous studies have developed new tools or used existing approaches to bias correct for other climate variables, such as the distribution-based scaling (DBS) approach for relative humidity and wind speed [71], the QM method for all climate variables [44], and the MCF method for solar radiation, relative humidity, and wind speed [54]. In this study, the QM-Gaussian method was also applied for other climate variables, including solar radiation, relative In this study, the QM-Gaussian method was used to correct the bias in maximum and minimum temperatures from GCMs and the results are very satisfactory. Although bias correction for temperature and precipitation is often made in other studies, bias correction for other climate variables is rarely conducted. Several previous studies have developed new tools or used existing approaches to bias correct for other climate variables, such as the distribution-based scaling (DBS) approach for relative humidity and wind speed [71], the QM method for all climate variables [44], and the MCF method for solar radiation, relative humidity, and wind speed [54]. In this study, the QM-Gaussian method was also applied for other climate variables, including solar radiation, relative humidity, and wind speed. But, the effect was not as good as that from temperature as the MCF method proposed by Haddeland et al., (2012) [54] is very simple and its bias correction effect had already been confirmed by other studies [72][73][74]. Therefore, the MCF method was adopted for the bias correction of other climate variables. As shown in Figure 7, the hydrological index RMSE was greatly improved after bias correction by the MCF method. Moreover, a Delta change approach was used to project the future PET in this study. Overall, the bias correction and Delta change method achieved the double effect of avoiding the possible impact of bias in climate variables from GCMs on PET estimation, which increased the accuracy of PET estimation in the future period.
As shown in Figure 8, the future PET values calculated by the FAO-PM, Makkink and Priestley-Taylor models from NorESM1_M (No.18) are much higher than other GCMs under the two RCPs. The common variable included in these three models but not included in temperature-based models is solar radiation. The future change of annual solar radiation from each GCM at the 12 stations are presented in Figure 15. In this figure, it can be observed that the variations of solar radiation from NorESM1_M are much bigger than those from the other 17 GCMs, which leads to the higher PET values. In addition, as shown in Figure 8, annual PET changes at the Linzhi and Longzi stations (No. 7 and No. 8) are smaller than other stations under NorESM1_M, which is owed to their smaller change of solar radiation.
Atmosphere 2019, 10, x FOR PEER REVIEW 21 of 25 QM-Gaussian method was also applied for other climate variables, including solar radiation, relative humidity, and wind speed. But, the effect was not as good as that from temperature as the MCF method proposed by Haddeland et al., (2012) [54] is very simple and its bias correction effect had already been confirmed by other studies [72][73][74]. Therefore, the MCF method was adopted for the bias correction of other climate variables. As shown in Figure 7, the hydrological index RMSE was greatly improved after bias correction by the MCF method. Moreover, a Delta change approach was used to project the future PET in this study. Overall, the bias correction and Delta change method achieved the double effect of avoiding the possible impact of bias in climate variables from GCMs on PET estimation, which increased the accuracy of PET estimation in the future period. As shown in Figure 8, the future PET values calculated by the FAO-PM, Makkink and Priestley-Taylor models from NorESM1_M (No.18) are much higher than other GCMs under the two RCPs. The common variable included in these three models but not included in temperature-based models is solar radiation. The future change of annual solar radiation from each GCM at the 12 stations are presented in Figure 15. In this figure, it can be observed that the variations of solar radiation from NorESM1_M are much bigger than those from the other 17 GCMs, which leads to the higher PET values. In addition, as shown in Figure 8, annual PET changes at the Linzhi and Longzi stations (No. 7 and No. 8) are smaller than other stations under NorESM1_M, which is owed to their smaller change of solar radiation. Data from 18 GCMs under the two RCPs were used to calculate PET in the future period and the uncertainty was also estimated in this study, which is beneficial for the analysis of the possible impact of climate change on PET and its contribution analysis. The results are helpful in understanding future PET changes and the selection of optimal PET models in the future period. The study of the appropriateness analysis of PET models under climate change is very important for many fields, including hydrologic modelling, agricultural water management, irrigation planning, crop yield simulation, drought monitoring and forecasting.

Conclusion
This study investigated the appropriateness of potential evapotranspiration models for climate change impact analysis, including one combination model (the FAO Penman-Monteith model) as a benchmark model, two temperature-based models (the Blaney-Criddle and the Hargreaves models) and two radiation-based models (the Makkink and the Priestley-Taylor models). The key findings of this study are summarized below.
Local calibration for the parameters of PET models is necessary, as performance evaluation indices (PBIAS and RMSE) are significantly improved after calibration. Data from 18 GCMs under the two RCPs were used to calculate PET in the future period and the uncertainty was also estimated in this study, which is beneficial for the analysis of the possible impact of climate change on PET and its contribution analysis. The results are helpful in understanding future PET changes and the selection of optimal PET models in the future period. The study of the appropriateness analysis of PET models under climate change is very important for many fields, including hydrologic modelling, agricultural water management, irrigation planning, crop yield simulation, drought monitoring and forecasting.

Conclusions
This study investigated the appropriateness of potential evapotranspiration models for climate change impact analysis, including one combination model (the FAO Penman-Monteith model) as a benchmark model, two temperature-based models (the Blaney-Criddle and the Hargreaves models) and two radiation-based models (the Makkink and the Priestley-Taylor models). The key findings of this study are summarized below.
Local calibration for the parameters of PET models is necessary, as performance evaluation indices (PBIAS and RMSE) are significantly improved after calibration.
In the baseline period, the PET estimations from the Blaney-Criddle and Makkink models are in better appropriateness with those obtained by the FAO-PM model. The Hargreaves and Priestley-Taylor models show a slightly poor performance.
The bias correction methods, i.e., the QM-Gaussian method for maximum and minimum temperatures and the MCF method for the other three climate variables, demonstrate their effectiveness by strongly decreasing RMSE after bias correction.
In the future period, monthly PET estimated by all models show positive trends. The changes of PET under RCP8.5 are much higher than RCP2.6. The radiation-based models (Priestley-Taylor and Makkink) show better appropriateness than the temperature-based models (Hargreaves and Blaney-Criddle) in the future, as the RMSE value of the former models is almost half of the latter models.