Periodic fluctuation of reference evapotranspiration during the past five decades: Does Evaporation Paradox really exist in China?

Evidence that the pan evaporation or reference evapotranspiration (ET0) as the indicator of atmospheric evaporation capability have decreased along with the continuous increase in temperature over the past decades (coined as “evaporation paradox”) has been reported worldwide. Here, we provide a nationwide investigation of spatiotemporal change of ET0 using meteorological data from 602 stations with the updated data (1961–2011). In addition, we explore the trigger mechanism by quantitative assessment on the contribution of climatic factors to ET0 change based on a differential equation method. In despite of different shift points regionally, our results suggest that the ET0 generally present decadal variations rather than monotonic response to climate change reported in previous studies. The significant decrease in net radiation dominate the decrease in ET0 before early 1990s in southern regions, while observed near-surface wind speed is the primary contributor to the variations of ET0 for the rest regions during the same periods. The enhancements of atmospheric evaporation capability after early 1990s are driven primarily by recent relative humidity limitation in China. From a continental scale view, as highly correlating with to Pacific Decadal Oscillation, the shift behaviors of ET0 is likely an episodic phenomenon of the ocean-atmosphere interaction in earth.

rising temperature during the recent two decades was identified by Liu et al. 23 with more recent observations. This phenomenon was further regionally verified in the hyper-arid region of Northwest of China by Li et al. 24 . Moreover, Wang et al. 25 suggested that there was actually a zigzag increasing-decreasing-increasing pattern of regional average ET 0 series with two joint points in 1973 and 1993 in the Tibetan Plateau. It is therefore worthwhile to re-examine the ET 0 patterns in larger area with different climatic and geomorphologic conditions based on the updated data. However, nationwide investigation on the changing patterns of ET 0 and its related trigger mechanism is not available.
The goal of this study is to examine whether the "evaporation paradox" is still true by investigating the response of ET 0 to changing climatic environment based on recently updated nationwide observation in China. Specifically, we begin by identifying the homogenously sensitive regions of ET 0 using a rotated empirical orthogonal function clustering method. Then we assess the spatial patterns of ET 0 trends and analyze the abrupt change in average series for individual homogenous regions. Meanwhile, special efforts are further made in characterizing the underlying aerodynamic and radiative driving mechanism for the ET 0 changes. Finally, the driving mechanisms under natural large atmosphere system are explored by identifying the correlation between changes in ET 0 and Pacific Decadal Oscillation.

Results
Identification of homogenous regions of ET 0 . To intensively characterize the spatial anomalies of ET 0 change nationwide, it is necessary to make a reasonable partition to the whole country. Traditionally, the whole China is always divided into eight climatic regions according to the latitude and longitude, and the climate regions are roughly coincide with the socioeconomic division 13,23 . However, there may be some uncertainty in sub-region selection based on visual inspection of the geographical distribution or administrative practices due to its vulnerability to subjective discretion 9 . This is especially true for the study of ET 0 changing patterns due to the fact that ET 0 involve in multifarious factors including geomorphologic features, topography factors and climate conditions. Therefore, the rotated empirical orthogonal function (REOF) method, a widely used clustering algorithm to objectively capture realistic spatial information, is employed in the current study to identify the homogenous regions regarding ET 0 in the whole China. The first 11 EOFs together explain 80.8% of the variance and the first 8 EOFs together can explain 73.7% of the variance. The fact that more than 70% of the total variability is captured in the first 8 EOFs indicates that the complexity of spatial patterns of ET 0 all over the country can mostly be explained by a small number of spatial structures. The percentages of variance of the new set of REOF modes generated by rotating the first eight loading vectors of the initial EOFs (Table S1) and correspondent isolines of the loading factor values ( Figure S1) suggest that Mainland China can be categorized into eight homogenous regions ( Fig. 1) by REOF analysis based on the annual ET 0 series from 602 stations for 1961-2011. Eight leading modes deduced sub-regions can roughly be described respectively as following: East Center China (EC), Southeast China (SE), Southwest China (SW), Tibet Plateau (TP), North Center China (NC), Northwest China (NW), North China Plain (NP), and Northeast China (NE). In term of the regional average series of these eight homogenous regions, investigation on temporal patterns and their attributions across the whole China are presented in following sections. Multidecadal changing patterns in ET 0 . We conduct trends and change point analysis for the area-averaged ET 0 in the resulted eight sub-regions by the means of the segmented regression model (shown in Fig. 2 and Table 1). The overall impression from the temporal patterns in ET 0 series is that shift trends can be evidently found in all eight sub-regions, although the turning years and the number of change points are regionally different. The decreasing-increasing (DI) patterns (with one joint point) are detected in EC, SE and SW regions with joint points happening at about 1992, 1992 and 1995, respectively. The zigzag increasing-decreasing-increasing (IDI) patterns (with two joint points) are depicted in TP (1973 and1993) and NW (1977 and1995) regions. With three joint points, approximating to periodic variation, more complex zigzag increasing-decreasing-increasing-decreasing (IDID) patterns are identified in NC (1973, 1991 and 2000) and NE (1978, 1993 and 2003) regions. It should be noted that although only two points of 1992 and 2002 are found in NP region, the overall changing pattern of decreasing-increasing-decreasing for ET 0 are similar to NC and NE regions, especially for the recent two decades. Interestingly, irrespective of complex shift trends reflected by different number of joint points and turning modes, multidecadal changing features in ET 0 are characterized by an evident spatial clustering partition. This clustering partition can roughly indicates that DI patterns mainly occur in southern region in China with a humid climate, whereas IDI mostly emerge in western Plateau regions with a cold and drought climate and IDID mainly appear in north and northeast regions with a cold climate. We also investigate the trends for the 602 stations during the segmented periods according to the corresponding regional change points. Taking EC region for example, as shown in Fig. 3, annual ET 0 at all stations decrease during 1961-1991 by 0-8 mm yr −1 (over 23% stations by 4-8 mm yr −1 ), of which over 80% stations are statistically significant at 95% confidence level (p < 0.05). As for as the second segmented period, from 1992-2011, more than 65% stations (over 30% stations by 4-12 mm yr −1 ) are dominated by increasing trends in ET 0 despite that only 35% stations are statistically significant at 95% confidence level (p < 0.05). Such extremely distinct spatial patterns in ET 0 during the adjacent segmented period are also identified in other homogeneous regions (shown in Figure S2-S8), giving us confidence in detecting multidecadal changes in ET 0 with shift trends in China.
Driving factors of changing ET 0 . The "evaporation paradox" is stated based on measured pan evaporation 3 , therefore, the correlation analyses of ET 0 with ET pan (data collected from183 stations with the same records of ET 0 ) is conducted before the contribution assessment of climatic factors. The close relationship between annual ET pan and ET 0 with correlation coefficient (R) of more than 0.9 suggests that ET pan is a good indicator of ET 0 across the whole China ( Figure S9). We identify the driving factors of changing ET 0 by quantifying the contribution of climatic factors to ET 0 change using partial derivatives (see Equation 10). For the eight homogeneous regions, the maximum relative error ρ(δ) representing the ratio of error δ to observed trend is found in NC region during the fourth segmented period of 2000-2011 with the value of −14.78%. The minimum ρ(δ) is only 1.70%, which occurs in NP region during the first segmented period of 1961-1991 (Table S2). The calculated ET 0 trends for all regions during all the segmented periods match very well with those detected from the computed ET 0 with P-M method (always regarded as observed ones) with quite satisfactory R 2 value of 0.99. The high accordance of calculated ET 0 trends with estimated ones reflected by high R 2 and low error indicates the differential equation method is applicable to quantify the contributions of climatic variables to changes in ET 0 . The contributions of climatic factors to the trends in ET 0 in categorized eight sub-regions during different periods are shown in Fig. 4. Taking NC regions with most complex shift trends in ET 0 as an example, during 1961-1972, the decreasing air temperature (T) causes a decrease of ET 0 at the rate of 0.42 mm yr −1 . Meanwhile, the increase of wind speed (U) and net radiation (R n ) as well as the decease of relative humidity (RH) lead to the increase of ET 0 at the rate of 2.53 mm yr −1 , 0.85 mm yr −1 and 1.18 mm yr −1 , respectively. Consequently, U is the dominant factor for the change in ET 0 in NC region during 1961-1972. Likewise, U, RH and R n are responsible for the decrease of ET 0 during 1973-1990, the increase of ET 0 during 1991-1999 and the decrease of ET 0 during 2000-2011 in NC region, respectively. Overall, although the change in ET 0 of the whole country is characterized by complicated spatial and temporal variability in the dominating contribution factors, some interesting phenomena can be summarized as below. Before the early 1990s, the decreases of ET 0 in EC, SE and SW regions, mostly located in South China, are mainly attributed to the decrease of R n . While for the rest regions including TP, NC, NW, NP and NE, mainly distributed in western and northern part of China, the changes of ET 0 before the early 1990s (all consist of increasing period before 1970s and subsequent decreasing period except NP region) mainly resulted from the changes of U, except the NE region during 1961-1977 with the dominating factor RH. From early 1990s to 2011, the ET 0 in EC, SE, SW, TP and NW regions present monotonous increase. While for the NC, NP and NE regions, changes in ET 0 are separated into two stages, i.e., increases from early 1990s to around 2000 and decreases after 2000. Although R n , T and U play respectively the most import role in the decrease of ET 0 at NC, NP and NE regions after around 2000, decreasing trends in RH is the most crucial factor for the increasing trends in ET 0 detected after the early 1990s in all the eight sub-regions.

The linkage between multidecadal patterns of ET 0 and the Pacific Decadal Oscillation (PDO).
In order to investigate the decadal characteristic of ET 0 in China and seek for the possible driving mechanisms, the regional mean anomalies of annual ET 0 , precipitation (P) and the climatic water availability (CWD, defined as P minus ET 0 ) as well as its relationship to the Pacific Decadal Oscillation index (PDO) are analyzed. Although as a whole, the regional ET 0 record shows a significant decrease trend of 0.298 mm yr −2 (p < 0.1) during the past 50-years period, the evolving process of continental mean ET 0 can be clearly divided into three stages. From 1961 to 1979, ET 0 present a slight but not significant decline (0.008 mm yr −2 ) with a high mean value. A prolonged significant decline starts in the year of 1980 and ends in the year of 1995, being contrary to the previously reported upward actual evapotranspiration in this period [26][27][28] . This may partly be explained by the complementary relationship 29 between ET 0 and actual evapotranspiration. After 1996, a pronounced recovery growth rate in ET 0 can be clearly revealed. An overall oscillation in regional mean ET 0 is found in China, which Scientific RepoRts | 6:39503 | DOI: 10.1038/srep39503 was not revealed by the previous studies, suggesting that the variation of ET 0 from 1961 to 2011 should be a phenomenon of earth's atmospheric circulation. PDO is considered as the important earth's climatic driving power for the continental hydro-meteorology abnormality such as shift of dry/wet in Asia, especially in China [30][31][32][33] . The PDO driven climate oscillations impart clear influences to the regional ET 0 in China as a whole, indicated by significant negative correlation between annual PDO index and annual mean ET 0 (r = − 0.384; p < 0.001). Years with high ET 0 mostly coincide with PDO cold phase although ET 0 responses to PDO may differ in different sub-regions in China. On the other hand, years with low ET 0 correspond to PDO warm phase conditions (Fig. 5). The regional mean P in China present large temporal variability with a negative but not significant trend of 0.257 mm yr −2 (p = 0.564) during 1961-2011, while the CWD present a significantly no trend although the slope is slightly positive with the value of 0.041 mm yr −2 (p = 0.938), suggesting that there is a weak tendency to reduce the deficit between atmospheric moisture requirement and available water supply for evapotranspiration. Despite of the strong relationship between regional ET 0 and PDO activities, the CWD is not significantly correlated with PDO (r = 0.164, p > 0.1) due to the weak relationship between P and PDO (r = 0.045, p > 0.1).

Discussion
The "evaporation paradox", reported in many regions of the world including China 15 , has drawn a considerable attention to explore the reason of the decreases in ET pan and/or ET 0 with increases in air temperature. However, with more recent observations or longer record, an increase in ET pan and/or ET 0 since 1980s or 1990s has been found in China or some regions of China 23,24,[34][35][36] . Specially, a zigzag increasing-decreasing-increasing pattern with two joint points in 1973 and 1993 has been identified in our previous study on the changes in ET 0 across the Tibetan Plateau 25 . However, the results from most previous studies on ET 0 change and driving mechanism are only based on subjective sub-regions selection. Meanwhile, the natural large atmosphere system triggering mechanisms for ET 0 change have not been investigated in previous studies. Consequently, a nationwide investigation in China is naturally performed in this study to detect the real multidecadal changing patterns of ET 0 and the correlation with PDO with extended data from 1961 to 2011. The temporal patterns of ET 0 in eight homogenous regions identified by REOF evidently reveal that there is no monotonous decreasing trend in ET 0 along with the continuously increasing temperature for the entire period (from 1961 to 2011), suggesting "paradox"  Table 1. Trend slope of reference evapotranspiration (ET 0 ) and meteorological variables in linear regression analysis, and the dominating factor result to ET 0 change by attribution analysis. Note: Trend slope of linear regression of ET 0 , T, RH, U, and R n are expressed in mm y −1 /y, °C/y, %/y, m s −1 /y, and W m −2 /y, respectively, and "*"represent significance at p = 0.10, "**"represent significance at p = 0.05.   24 . When other sub-regions are considered, coarse periodic behaviors of ET 0 change with two/or three turning point seems to be reasonable as compared with total decreasing-increasing patterns of evaporation in entire China reported by Cong et al. 38 and Liu et al. 23 . The combined effects of climatic variables to ET 0 changes are revealed in all the sub-regions with the sensitive analysis, and are effective with a good agreement between the observed and calculated ET 0 trends (see Table S2). Although relative humidity is the most sensitive variable (followed by net radiation, wind speed and air temperature) for all the period as well as all the sub-regions except SW and TP region, there are different dominating factors for ET 0 change because the contributions of climate factors depend on the combined effect between sensitivity and trends of a variable itself 9 . Decline of net radiation only plays the important role for the ET 0 decrease before early 1990s rather than the entire term in EC, SE and SW regions located in south China, highlighting the potential implications of "global dimming" in China over a roughly 1960-1990 period and a recovery (coined with "global brightening") thereafter 39,40 , which is consistent with previous studies in China 13,20,23,38,41 , the United States 3 , the northeast of India 42 , and Greece 43 . There is consensus among many researchers that cloud coverage and aerosol, which are not completely independent variables with interaction in various ways 44 , are considered as the most likely candidates for the explanation of global dimming and brightening 39,[45][46][47][48] . Particularly strong evidence for aerosol effects on surface solar radiation with increasing air pollution in China were noted in various studies [49][50][51][52] . As far as the other regions with larger inter-annual variability in ET 0 are concerned, variation of wind speed generally plays more decisive role in the change of ET 0 from 1961 to early 1990s, in line with the evidences reported in many places and summarized by McVicar et al. 11 . Although the dynamic mechanism of recent slow-down in near-surface global winds is complicated, it can generally be attributed to the increases of terrestrial surface roughness 53,54 or large scale atmospheric circulations [55][56][57][58] . Focusing on the periods after 1990s, our preliminary analyses have shown that dominating factors for the increases of ET 0 transited from net radiation and wind speed to relative humidity for all the sub-regions. This coincides with the recent regional study in the arid regions of northwest china by Li et al. 24 and national wide study in China by Cong et al. 38 . A global evidence of reduction in surface atmospheric humidity revealed in this study in national scale was presented by Simmons et al. 59 . They inferred the recent reduction in relative humidity over land may be due to limited moisture supply from the oceans. For the purpose of simplification, the systemic errors of contribution analysis reflected by the differences between observed ET 0 and complied trends may be due to the interactions among climatic factors (which were not considered here but may exist in practice) and ignorance of other factors such as aerosol and dust. From global perspective, there is an indication for a shift that occurred in the late 1970s with minor change occurred around 1990 and 2000 for PDO 60,61 . Having the strong correlation with PDO, the decadal variations of ET 0 in China with two turning points during 1961-2011 are likely dominated by the episodic dynamics of climatic system, indicating that the recent ET 0 trends in China reflect PDO conditions and are not the consequence of a persistent reorganization of the hydro-meteorological cycle. The similar conclusions have also been obtained from the global analysis for continental evaporation 27,28 . Consequently, the widely reported "evaporation paradox", a counterintuitive behavior, representing that evapotranspiration decrease despite the increase of air temperature, seem to be more reasonably explained as part of a climate oscillation, at least in China 33 . The emerging picture of enhanced of ET 0 in southern China highlights the possible threat posed by acceleration of the terrestrial hydrological cycle to water resources management and food security in an enhanced-greenhouse-affected climate.

Materials and Methods
Data. Daily measurements of air temperatures (minimum, maximum and average) at 2 m height, relative humidity, wind speed at 10 m height, and sunshine duration obtained from the National Climatic Centre (NCC) for 602 ground-based stations in China during the period 1961-2011, provided by the National Meteorological Information Centre of China (NMIC) of the China Meteorological Administration (CMA), are used for estimating ET 0 by Penman-Montieth method. Wherein, wind speed was adjusted to 2 m height by virtue of wind profile relationship proposed by Allen et al. 62 . Net radiation at 51 stations was available. Moreover, precipitation data from these meteorological stations were also collected.

Food and Agriculture Organization (FAO) Penman-Montieth (P-M) Method.
From the reference surface defined as "a hypothetical reference grass with an assumed crop height of 0.12 m, a fixed surface resistance of 70s m −1 and an albedo of 0.23", the daily reference evapotranspiration (ET 0 ) is estimated by the Penman-Montieth (P-M) combination method, which is recommended by the FAO as a standard to calculate ET 0 wherever the required input data are available and is proven to have good performance for various climatic conditions worldwide 9 . The P-M equation can be expressed as (Allen et al.) 62 : where ET 0 is the reference evapotranspiration (mmd −1 ), R n is the net radiation at the crop surface (MJ m −2 d −1 ), G is the soil heat flux density (MJ m −2 d −1 ), T is the mean daily air temperature (°C), u 2 is the daily average wind speed at 2 m above ground level (ms −1 ), e s is the saturation vapor pressure (kPa), e a is the actual vapor pressure (kPa), e s − e a is the saturation vapor pressure deficit (kPa), Δ is the slope of the saturated vapor pressure in relation to air temperature (kPa °C −1 ) and γ is the psychrometric constant (kPa °C −1 ). As the difference between incoming net shortwave radiation (R ns ) and outgoing net longwave radiation (R nl ), R n can be expressed as: n n s n l R ns is derived from the balance between incoming and reflected solar radiation and is given by: where a is the albedo or canopy reflection coefficient, which is 0.23 for the hypothetical grass reference crop. R s can be estimated from sunshine duration (or hours of sunshine) with the help of the Angstrom formula where R s is the solar or shortwave radiation (MJ m −2 d −1 ), n is the actual duration of sunshine (h), N is the maximum possible duration of sunshine or daylight hours (h), ( n N is thus the relative sunshine duration), R a is the extraterrestrial radiation (MJ m −2 d −1 ), a s and b s are the Angstrom coefficients. We calibrate the Angstrom coefficient using the observation of n and R s from 51 stations widespread the whole China. Meanwhile, Ra and N are computed based on date and latitude according to the equations provided by the FAO P-M method. Overall, the high R 2 values (see Table S3) between observed and estimated R s indicate that the Angstrom model is suitable for daily global radiation estimation in China. For stations with no observation of solar radiation but sunshine duration, a s and b S are estimated by Kringing interpolation method.
R nl is based on water vapor, clouds, carbon dioxide and dust are absorbers and emitters of longwave radiation, and can be estimated by where T max,K and T min,K are respectively maximum and minimum absolute temperatures during the 24-hour period, K = °C + 273. 16. The soil heat flux, G, is the energy that is utilized in heating the soil. G is positive when the soil is warming and negative when the soil is cooling. The soil heat flux is small compared to Rn and may often be ignored in daily evapotranspiration estimating.
Trend analysis and breakpoint identification. A simple linear regression method was used to calculate the slope of linear least squares regression line fit to the inter-annual variation of ET 0 and other climatic variables inputted in the P-M equation. The significance of changes of these variables across China was mapped and assessed with the help of the two-tailed significance tests. The segmented regression with constraints method developed by Shao and Campbell 63 , which can detect both shift trends and step changes simultaneously without knowing the number of trend segments and change points and their location over time 64 , is employed to investigate the possible turning points of the regional average ET 0 series across China. The method was summarily described as follow 63-65 : Assuming that x i means the observed variable at time t i (i = 1, 2, … , N) t i , the regression can be written as where the errors {ε i } (i = 1, 2, … . N) are assumed to be independent and identically distributed as N(0, σ 2 ). For L abrupt change points r 1 < r 2 < ···r L , resulting in L + 1 segments, the model can be written as l l l 1 with r 0 = 0, r L+1 = ∞ and u l (t) are the deterministic part quantifying the trends. Suppose that there are j l join points {J l., k } (k = 1, 2, … . j l ) in the ith segment between the break points r l−1 and r l , u l (t) can be defined as: A modified Akaike's information criterion (AICc), derived by Hurvich and Tasi 66 with c standing for the second order correction 67 , is employed to select an optimal numbers of break and join points and their locations, and the optimal model can be obtained by minimizing AICc.

The Rotated Empirical Orthogonal Function (REOF) method.
To characterize in detail the spatial variability of ET 0 at nationwide scale and identify the homogenous regions, we applied REOF to analyze the most dominant spatial patterns. The aim of the EOF method is to find a relatively small number of independent variables conveying as much of the original information as possible without redundancy by decomposing a multivariate data set into uncorrelated linear combination of separate function of the original variables. However, the physical interpretability of the obtained patterns is sometimes a matter of controversy because of orthogonality in both space and time 68 . This limitation has brought about the development of the rotated empirical orthogonal function (REOF) 69 , which can cluster within each mode a small number of high valued variables and a large number of near-zero value variables through a rotation to simple structure. In a comparison study by Kim and Wu 70 , REOF was found to be better in dividing climatic patterns. In this paper, with maximizing the variance of the squared correlation between each rotated principal component (RPCs) and each variable, the varimax REOF method was chosen to give the simplest pattern description while explaining the maximum amount of variance. Assuming that y and x i are the time series variables, dy dt and dx dt i should be the long-term trend for y and x i , respectively. The term of ′ f i dx dt i , the product of long-term change in x i and the partial derivative, can then be considered as the contribution of change in x i to the long-term variation of y. Thus, following P-M equation for ET 0 estimation, contributions of changes in key climatic factors to ET 0 variation can be approximately as 55   where S(x i ) is the sensitivity coefficient of ET 0 related to x i , the i th variable. Being first employed by McCuen 71 , the sensitivity coefficient has been widely used in evaluating the climatic response of evapotranspiration, especially during recent two decades [73][74][75][76][77] .