Changes in reference evapotranspiration over China during 1960–2012: Attributions and relationships with atmospheric circulation

This study investigates reference evapotranspiration (ET0) trends in China from 1960 to 2012 based on the Penman–Monteith equation and gridded meteorological measurements. Under the combined impacts of factors influencing ET0 (i.e., net radiation [RN], mean temperature [TAVE], vapour pressure deficit [VPD], and wind speed [WND]), both seasonal and annual ET0 for the whole China and more than half of the grids decreased over the past 53 years. The attribution analyses suggest that for the whole China, the WND is responsible for annual and seasonal ET0 decreases (excluding summer, where RN is responsible). Across China, the annual cause of WND with the largest spatial extent (43.1% of grids) mainly derives from north of the Changjiang River Basin (CJRB), whereas VPD (RN) as a cause is dispersedly distributed (within and to the south of the CJRB). In summer, RN is dominant in more than half of the grids, but the dominance of VPD and WND accounts for approximately 90% of grids during the remaining seasons. Finally, the correlation coefficients between ET0 and the Atlantic Oscillation (AO), North AO, Indian Ocean Dipole (IOD), Pacific Decadal Oscillation (PDO), and El Niño Southern Oscillation (ENSO) indices with different lead times are calculated. For the whole China, annual and seasonal ET0 always significantly correlate with these indices (excluding the IOD) but with varied lead times. Additionally, near half of the grids show significant and maximum (i.e., the largest one between ET0 and a certain index with a lead time of 0–3 seasons) correlation coefficients of ET0 with PDO in spring and summer, ENSO in autumn, and AO in winter. This study is not only significant for understanding ET0 changes, but it also provides preliminary and fundamental reference information for ET0 prediction.

droughts, and extreme precipitation) have become a focus of recent climate change studies (Milly, Wetherald, Dunne, & Delworth, 2002;Mishra, Ganguly, Nijssen, & Lettenmaier, 2015;Sheffield & Wood, 2008;Trenberth, 2010;Zhang, Cheng, Zhang, & Liu, 2017). There is a concern that with climate change, droughts may increase in intensity, frequency, area, and duration (Dai, 2013), leading to great impacts on human and natural systems and providing an intense reminder of the potential consequences of droughts. Recently, based on the close relationship between drought and precipitation (Heim, 2002), numerous studies have completely explained drought mechanisms across the globe (Buttafuoco, Caloiero, & Coscarelli, 2015;Schubert et al., 2016;Yan, Zhang, Zhou, & Han, 2017). However, it should be noted that the impacts of evapotranspiration (ET) or reference ET (ET 0 ) have been given less attention, even though they are closely associated with changes in climatic variables and are integrated indicators of climate change. For example, Hu and Willson (2000) and Sun et al. (2016),  have reported that the impact of ET 0 on drought could be equivalent to or even higher than that from precipitation. Consequently, understanding the spatial and temporal variations in ET 0 is a vital component for comprehensively understanding hydrometeorological changes and is helpful for better monitoring and predicting drought (Beguería, Vicente-Serrano, Reig, & Latorre, 2014;Li, He, Quan, Liao, & Bai, 2015;McEvoy, Huntington, Mejia, & Hobbins, 2016;Sun et al., 2016;; Vicente-Serrano, Beguería, & Lópezmoreno, 2010;Zhang et al., 2018).
Climate change is intensifying and mainly characterized by significant increases in temperature (IPCC, 2014). However, ET 0 has not increased with global warming as expected. In contrast, it has shown a decreasing trend in most regions of the world (Roderick & Farquhar, 2002), which is known as the evaporation paradox (Brutsaert & Parlange, 1998). To explain this issue, a number of studies have been performed by quantifying the impacts of meteorological factors on ET 0 , and they concluded that changes in wind speed and radiation (or sunshine duration) are the main causes of decreased ET 0 in some regions (Chen, Liu, & Thomas, 2006;Jhajharia, Dinpashoh, Kahya, Singh, & Fakheri-Fard, 2012;Peterson, Golubev, & Groisman, 1995;Roderick, Rotstayn, Farquhar, & Hobbins, 2007;Stanhill & Möller, 2008), whereas water vapour indicators (e.g., relative humidity and vapour pressure deficit) and temperature exert dominant impacts over other regions (Liu, Zheng, Liu, & Cao, 2009;Zhang, Liu, & Hong, 2013).
Recently, to overcome this issue (at least partially) and enhance the confidence level of the conclusions,  developed a new separation method to attribute ET 0 changes over Southwest China and stated that this method has a better performance (e.g., accuracy) compared with other methods. As a result, this provides an efficient tool to more accurately and comprehensively understand the underlying mechanisms of changed ET 0 across China. Furthermore, by considering the importance of ET 0 for the estimation of crop water demand and drought forecasting systems (McEvoy et al., 2016), obtaining detailed information about ET 0 is critical before developing countermeasures for policy makers. Large-scale atmospheric or coupled atmosphere-ocean modes of climate variability (usually referred to as teleconnections due to their long-distance range of influence) can provide a useful framework for linking ET 0 to climate fluctuations via influencing climate variables (e.g., wind speed, solar radiation, temperature, and relative humidity; Chen, Li, & Pryor, 2013;Li & Chen, 2014;Yu, Zhong, Bian, & Heilman, 2015;Zhu et al., 2017). Given more knowledge of atmospheric system processes and the improved capability of future teleconnection predictions (Alexander, Matrosova, Penland, Scott, & Chang, 2008;Derome, Lin, & Brunet, 2005;Yang et al., 2015;Zheng, Fang, Zhu, Yu, & Li, 2016), estimating future ET 0 has become more possible using historical and predicted teleconnection patterns. However, few studies (Meza, 2005;Sabziparvar, Mirmasoudi, Tabari, Nazemosadat, & Maryanaji, 2011;Tabari, Talaee, Some'e, & Willems, 2014) Tabari et al. (2014) found that winter ET 0 over Iran was negatively correlated with the North AO index at almost all of the chosen stations. In summary, comprehensive analyses of the linkages between ET 0 and teleconnection patterns can provide important information for ET 0 prediction.
In short, in-depth discussions about the causes of ET 0 changes on multiple timescales (e.g., seasonal and annual), and the linkages between ET 0 and teleconnection patterns can provide important informational support for accurately estimating and predicting agricultural water requirements and reasonably making irrigation decisions, as well as facilitating a deep understanding of regional dry/wet (drought) conditions (Meza, 2005;Saadi et al., 2015;Shi et al., 2014). Therefore, the main objectives of this study are as follows: (a) to analyse seasonal and annual ET 0 changes across China from 1960 to 2012, (b) to quantify the respective impact of each meteorological factor on ET 0 changes and identify the dominant factors, and (c) to explore the relationships between teleconnection patterns and ET 0 .
2 | DATA AND METHODOLOGY

| Data
To compute ET 0 , monthly maximum and minimum temperature (°C), relative humidity (%), sunshine duration (hr), and wind speed at a 10-metre height (m/s) are measured from 1960 to 2012 at 1672 weather sites, which are collected from the China Meteorological Administration. These sites and the 10 basins across China are shown in Figure 1. Before using, two data quality issues within this dataset (i.e., inhomogeneous and missing values) should be noted; therefore, we have followed the methods of  to solve them. For simplicity, the detailed procedure is not presented here but can be found in . Considering that a nonuniform density distribution of the sites has potential impacts on the results, we follow J.  and process 1,209 sites into 548 grid boxes of 1°× 1°(longitude by latitude). Each grid box should contain at least one site with continuous meteorological measurements. If any grid includes more than one site, the average value of each meteorological variable at these sites is calculated to represent the final value of that grid. Notably, spring, summer, autumn, and winter are specified as March-May, June-August, September-November, and December-February, respectively.
China is located in a typical monsoon region (i.e., the East Asian monsoon region), which is significantly influenced by ENSO (Ding et al., 2014;He & Wang, 2013). The PDO can also regulate China's climate (Feng, Wang, & Chen, 2014;Wang, Chen, & Huang, 2008). In addition, the East Asian monsoon has been influenced by the IOD and AO (Chen, Feng, & Wu, 2013;Ding, Ha, & Li, 2010;He, Gao, Li, Wang, & He, 2017;Yuan, Yang, Zhou, & Li, 2008 We compute ET 0 using a modified Penman-Monteith equation (Allen, Pereira, Raes, & Smith, 1998), which has been recommended as a standard tool for climatic data by the Food and Agriculture Organization of the United Nations and can be written as follows: (1) where RN (MJ/(m 2 d)) is the difference between net incoming shortwave and net outgoing longwave radiation; G (MJ/(m 2 d)) denotes the soil heat flux, which can be ignored at monthly or longer timescales; γ (kPa/°C) and Δ (kPa/°C) represent the psychrometric constant and slope vapour pressure curve, respectively; TAVE (°C) represents the mean temperature and can be estimated as the mean of the minimum and maximum temperature; WND (m/s) is the wind speed at a 2-metre height, which is converted from the wind speed at a 10-metre height; and VPD (kPa) denotes the vapour pressure deficit, which can be expressed as the difference between saturation and actual vapour pressure. All of the computations are performed on a monthly scale. Detailed equations for Δ, RN (which can also be found in Supporting Information), γ, TAVE, WND, and VPD can be found in Allen et al. (1998).

| Method for attributing ET 0 changes
To quantify each the individual contribution of each factor to the ET 0 trend, five experiments are designed, namely, one control experiment Then, to distinguish the contribution from each factor alone, we can utilize the following equation: where A C i represents the contribution of the i-factor to ET 0 change, and T Sim_CTR and T Sim_i denote the ET 0 trend for Sim_CTR and Sim_i, respectively. Here, this method is named as approach A.
Due to the interactions among these driving factors, there exists the potential to introduce some uncertainties into the attribution of ET 0 change. To eliminate the complex effects, Sun et al. (2014) and  proposed a new separation method (i.e., approach B), which is shown as follows: where ∑ n k≠i B C k represents the accumulative contributions of all of the driving factors (except for the i-factor) to the ET 0 trend; n is equal to 4 here and indicates the number of sensitivity experiments, and T Sim_i represents the ET 0 trend via Sim_i. As a result, the contribution of each factor alone can be obtained by solving Equation 3: For greater robustness and accuracy of the results, the performances of both approaches A and B are evaluated, which are shown in Table 2 and Figure S1. For the annual ET 0 trend, the root mean square error (RMSE; mean relative error, MRE) is 0.057 mm/year (51.62%) and 0.019 mm/year (17.21%) for approaches A and B, respectively. Table 2, all seasonal RMSEs (<0.01 mm/year; MREs < 18%) via approach B are always smaller than those via approach A. In addition, Figure S1 displays the scatterplots of the cumulative contributions (∑ 4 k¼1 L C i , where L is A or B) and the Sim_CTR ET 0 trend at annual and seasonal scales. It is evident that relative to approach A, the results via approach B are closer to the 1:1 line ( Figure S1). Therefore, approach B (with a higher accuracy and efficiency) is chosen to quantify the impacts of climate change on the ET 0 trend in this study.

| Relationships between ET 0 and teleconnection indices
In this study, correlation coefficients are used to measure the strength  the SHRB and LRB in spring). In addition, it is worth noting that over the SHRB and SWRB, ET 0 increases from summer to winter, particularly for ET 0 during summer and autumn over the SWRB, which has larger (approximately 0.3 mm/year) and more significant trends.
As shown in Figure 3a and Table 3 Table 3), grids with increased (decreased) ET 0 account for FIGURE 2 Seasonal and annual ET 0 trends from regional time series of each basin. Asterisk denotes that the trend is significant with P < 0.05 FIGURE 3 Spatial distributions of seasonal and annual ET 0 trends. Asterisk denotes that this trend is significant with P < 0.05 approximately 50% of the total, with significant changes generally in the mideastern CJRB and northern SERB (eastern SHRB, LRB, and HRB, and north-western HaRB). During summer (Figure 3c and  (Figure 3d and  (Figure 3e) is consistent with that in summer ( Figure 3c), which is characterized by a decrease in most of the grids (61.6% in

| Changes in major meteorological factors
Changes in ET 0 are jointly affected by various meteorological factors, such as TAVE, VPD, WND, and RN; thus, it is necessary to know their variations to understand ET 0 changes. Figure 4 presents the seasonal and annual trends of these four factors in China and its 10 major basins. For the China regional average (the right panel of Figure

FIGURE 4
Seasonal and annual trends of ET 0 driving factors from regional time series of each basin and the whole China. Asterisk denotes that the trend is significant with P < 0.05 Over China and each basin (Figure 4b), the regional average seasonal  Figure 4d). In addition, we calculate the annual and seasonal changes in these driving factors, which are illustrated in Figures S4-S7. Overall, the annual and seasonal TAVE and VPD (RN and WND) increase (decrease) for most (>70%) of the grids (Table S1), whereas the changing rates have obvious spatial differences.

| Causes of ET 0 change
Using the separation method proposed by Sun et al. (2014) (Table 4). On an annual scale (Figure 7a and

| Relationships between ET 0 and various teleconnection indices
Annual correlation coefficients between ET 0 and the four selected teleconnection indices are first calculated for the whole China and all of the basins, which are presented in   Table   S2). Significant AR et0,PDO values (Figure 8(a3) and Table S2) are found in 28.6% of grids, which basically cover the NWRB, middle CJRB, and eastern PRB. In general, AR et0 , ENSO is smaller across China, and only several grids have significant correlations (Figure 8(a4)). As shown in Figure 8(b1-b4), the spatial extent of significant AR et0,y_1 obviously decreases compared with that of significant AR et0,y (except for that influenced by ENSO; Table S2). From the perspective of spatial distributions, AR et0,AO_1 (AR et0,PDO_1 ) has a similar spatial pattern with that of AR et0,AO (AR et0,PDO ), which is generally characterized by significant   Figure 10 shows the spatial patterns of the lead season (s) for each index, which shows the MSR with seasonal ET 0 . In spring (Figure 10(a1) and  (Table   S3). It is evident that during summer (Figure 10(b1-b4) and  Note. Bold number denotes that the correlation coefficient is significant with P < 0.05.

FIGURE 8
Spatial distributions of annual correlation coefficients between ET 0 and four teleconnection indices. Asterisk denotes that the correlation coefficient is significant with P < 0.05 FIGURE 10 Spatial distributions of lead season (s) of each index, which corresponds to the MSR. Asterisk denotes that this MSR is significant with P < 0.05

FIGURE 9
Maximums of seasonal correlation coefficients between ET 0 and a certain index with a lead time of 0-3 seasons for each basin and China. Bold number suggests that the correlation coefficient is significant with P < 0.05 HaRB, north-central CJRB, and middle PRB), but for the AO and PDO, these grids appear in regions north of the CJRB (Figure 10(c1-c4)).
Regarding significant winter MSRs (Figure 10(d1-d4) and Table S3), grids influenced by the IOD, PDO, and ENSO are limited (<18% of grids) and basically occur in northern SHRB, eastern LRB, and the mouth of the CJRB, SHRB, NWRB, and southern SWRB, and SERB, eastern PRB, and southern SWRB, respectively. However, 49% of grids with significant MSRs are influenced by the AO and widely distributed across China, including the northern SHRB, YRB, PRB, SERB, southern SWRB, and middle CJRB. In addition, based on comparisons of the spatial distributions and grid percentages of significant MSRs ( Figure 10 and Table S3), we can conclude that the PDO is more significant for predicting spring and summer ET 0 for approximately half of the grids, but ENSO and the AO influence autumn and winter ET 0 estimates, respectively.

| Comparison of results
In the present study, a decrease in WND (RN) is concluded to be a major determinant factor for annual ET 0 changes for the whole China

| Causes for trends in meteorological factors
Similar to those in other parts of the globe (Roderick & Farquhar, 2002;Stanhill & Cohen, 2001;Vautard, Cattiaux, Yiou, Thépaut, & Ciais, 2010), both seasonal and annual WND and RN are found to decrease from 1960 to 2012 for the whole China and most of the 548 grids; furthermore, the findings generally coincide with the previous results in China and its subregions (e.g., Lin, Yang, Qin, & Fu, 2013;Wang, Yang, Han, Wang, & Zhang, 2013;Xu et al., 2006;Yang, Zhao, Hu, & Zhou, 2009). To explain the decreased WND, many scholars have conducted various studies from the perspective of large-scale atmospheric circulation (e.g., Jiang, Luo, Zhao, & Tao, 2010;Xu et al., 2006). For example, Jiang et al. (2010) and Xu et al. (2006) pointed out that decreases in temperature vary between the Asian continent and the low-latitude region of the ocean, which are responsible for decreased WND over China. In addition, human activities, such as land use/cover change (e.g., urbanization) and air pollution, can also exert impacts on WND (Ren et al., 2008;Wang, Feng, & Gao, 2014;Xu et al., 2006). Regarding the possible causes of RN reductions, previous studies have suggested that cloud coverage and air pollution and aerosol loading from rapid urbanization are potential reasons (Qian, Kaiser, Leung, & Xu, 2006;Xia, Wang, Chen, & Liang, 2006;Yang et al., 2009). Zhu et al. (2017) found that large-scale atmospheric or coupled atmosphere-ocean modes are strongly associated with cloud cover across the world, which leads to RN variations. Moreover, with a weaken WND in China, air pollution and aerosol loading cannot be rapidly dispersed, which results in the further decrease in RN (Yang et al., 2009).
Under the background of global warming possibly caused by anthropogenic activities (e.g., greenhouse gases emissions; IPCC, 2014;Xu, Gao, Shi, & Zhou, 2015), annual and seasonal TAVE basically respond to significant increases at national and grid scales in our study. Notably, there are evident differences in warming rates across China, which are closely associated with the regional characteristics of the land surface's energy balance (Dong, Xi, & Minnis, 2006;Stanhill & Ahiman, 2014;Wild, Grieser, & Schär, 2008;Wild, Ohmura, & Makowski, 2007), especially for the longwave radiation budget. Notably, temperature fluctuations are also linked to climate variability signals (Chen, Chen, Bai, & Xu, 2016;Guo, Zhao, & Dong, 2015). For example, Guo et al. (2015) noted that the AO plays an important role in variations in winter surface air temperature over China. Chen et al. (2016) found that the PDO can exert significant impacts on interdecadal temperature fluctuations over Northwest China, whereas the IOD plays a more important role in annual warming. For an overwhelming majority of grids, there are consistent increases in annual and seasonal VPD, which is consistent with the conclusions from Yang and Yang (2012). It is well known that the VPD is jointly controlled by TAVE and relative humidity (RH); furthermore, the VPD experiences increases (decreases) with increased TAVE (RH) when RH (TAVE) is constant and vice versa. To address this issue, we have analysed RH changes across China and found that annual and seasonal RH generally decreases across China (not shown here), which is possibly due to less moisture from oceans (Li & Chen, 2014;Simmons, Willett, Jones, Thorne, & Dee, 2010). Thus, we conclude that increased VPD for most of China can be attributed to increased TAVE and decreased RH.
Additionally, based on the discussions above and previous studies (e.g., Chen et al., 2016;Guo et al., 2015), we can basically obtain the underlying mechanisms of teleconnection impacts on ET 0 by causing variations in climate variables (e.g., TAVE, RN, relative RH, and WND). This is why teleconnection indices are significantly correlated with ET 0 at different timescales in some regions (Figures 8-10).  Reid et al. (2003) found that due to the impacts of influential elements (e.g., water supply) associated with the physiological functions of plants, the stomatal characteristics of plants do not change with increasing CO 2 . Li, Kang, and Zhang (2004) and Peng, Dan, and Dong (2014) reported that increases in CO 2 are conducive to plant growth and leaf area increases, which consequently intensify ET. In a word, elevated CO 2 impacts on ET are still questionable to date. Land-surface albedo (set as a constant of 0.23 here) change is also a source for uncertainties in this study. This parameter determines the ability of land surfaces to absorb solar radiation, which can affect the ET process by altering land-surface radiative forcing (Twine, Kucharik, & Foley, 2009). As we know, albedo differs among various types of land use/cover (Barnes & Roy, 2008); therefore, land use/land cover changes likely change albedo (Lee et al., 2011;Wu, Zha, & Zhao, 2015). In recent years, due to rapid social and economic developments in China, major adjustments and structural changes have occurred in the spatial patterns of land use/cover (Liu et al., 2010). It has been reported that recent land use/land cover changes in China have resulted in a general decrease in albedo Zhai, Liu, Liu, Zhao, & Huang, 2014) but with obvious regional differences. Taking the Beijing-Tianjin-Tangshan area as an example, the conversion of agricultural land to urban areas decreases albedo, which results in more energy for ET (Zhai et al., 2014). However, less energy in the Sanjiang Plain is partitioned for ET due to increased albedo, which is mainly associated with the replacement of forests by farmlands (Zhang, Wang, Fang, Ye, & Zheng, 2012). Land-surface resistance can control the ET process by influencing the flow resistance of water vapour during crop transpiration and surface evaporation (Allen et al., 1998;Boegh, Soegaard, & Thomsen, 2002). Shuttleworth and Gurney (1990) pointed out that land-surface resistance is associated with the plant type, climate conditions, canopy structure and plant-soil water status, implying that land-surface resistance varies with space and time. However, the impacts of these factors on resistance are not included in the Penman-Monteith model, as this parameter is set to 70 s/m based on the most suitable condition without water stress (Gharsallah, Facchi, & Gandolfi, 2013). As a result, a fixed land-surface resistance tends to produce uncertainties. In spite of the used separation method with a better performance (i.e., lower MRE and RMSE), biases between the cumulative contributions from all driving factors and the Sim_CTR trends still exist ( Figure S1); therefore, we should keep in mind the uncertainties from the limitations of this method. First, the key concept of this separation method is that ET 0 changes are linearly and jointly induced by all factors; however, each factor's impacts are To explore the linkage between ET 0 and teleconnection patterns, we have calculated linear correlation coefficients between ET 0 and four indices with different lead times. It is beyond a doubt that the analyses provide preliminary and necessary information (i.e., which index has a significant correlation with ET 0 and where) for predicting ET 0 . Nevertheless, it is a fact that the relationships between ET 0 and teleconnection patterns are complex and nonlinear; therefore, using linear correlation analyses is not sufficient for understanding the linkage between ET 0 and these indices. In past decades, general atmospheric circulation models or earth system models have made great progress and provided ideal and efficient tools to comprehensively study the underlying mechanisms of various hydrometeorological

| Uncertainties
issues. However, it should be noted that our major aims are to attribute ET 0 changes and quickly identify potential predictors for estimating ET 0 , and there are still some difficulties when diagnosing the linkages between ET 0 and teleconnection patterns with models (i.e., how to reasonably design numerical experiments due to the higher complexity of the atmosphere system). Therefore, applying models for investigating the relationships between ET 0 and these indices may be beyond the scope of this study (at least now). Regardless, we would like to conduct in-depth studies about this issue using models in the future, which will enhance the confidence level of our conclusions.

| CONCLUSIONS
Comprehensive To sum up, this detailed information about the attributions of ET 0 changes will be useful to further understand dry/wet conditions (e.g., ongoing and intensifying droughts over some parts of China) and guide policy makers for the management of water resources and agricultural production activities.
Lastly, annual and seasonal correlations between ET 0 and the four teleconnection indices (i.e., IOD, AO, PDO, and ENSO) with different lead times are estimated. In general, there are always significant correlations between ET 0 and each index in some regions, but the locations and spatial extents vary considerably among different indices at given lead times. Overall, the AO and PDO have the potential to be predictors for forecasting annual ET 0 at different spatial scales (i.e., the whole China, basin, and grid scales). For each season, the PDO has greater significance for predicting spring and summer ET 0 , but ENSO and AO are more significant for predicting autumn and winter ET 0 , respectively. Undoubtedly, these correlation analyses will provide fundamental and important reference information for predicting ET 0 with the statistical forecast method, which has a potential application in climate forecasting systems.