Trend, Seasonal, and Irregular Variations in Regional Actual Evapotranspiration Over China: A Multi-Dataset Analysis

Actual evapotranspiration (AE) is a crucial processes in terrestrial ecosystems. Global warming is expected to increase AE; however, various AE estimation methods or models give inconsistent trends. This study analyzed AE variability in China during 1982–2015 based on the Budyko framework (AE_Budyko), a complementary-relationship-based product (AE_CR), and the weighted average of six reanalyses (AE_WAR). Because the response of AE to driving factors and the performances of AE datasets are both scale-dependent, China has been categorized into six distinct climatic areas. From a regional perspective, the X-12-ARIMA method was used to decompose monthly AE into the trend, seasonal, and irregular components. We examined the main characteristics of these components and the relationships of climate factors with AE. The results indicate that the trend component of AE increased from the mid-1990s to the early 2000s and more recently in the hyper-arid and arid areas. Increasing AE was observed from 1982 to the early 1990s in the semi-arid and dry sub-humid areas. AE increased significantly and had substantial interannual variability for the entire period in the sub-humid and humid areas. Increased precipitation and water supply from terrestrial water storage contributed significantly to increasing AE in the drylands. The simultaneous occurrence of increasing precipitation and wet-day frequency caused increasing AE in the dry sub-humid area. Increased AE could be explained by the increased energy supply and precipitation in the sub-humid and humid areas. Precipitation had the strongest influence on the irregular component of AE in drylands. AE and potential evapotranspiration had a strong positive correlation in the sub-humid and humid areas. Regarding data availability, a discrepancy existed in the trend component of AE_CR because soil moisture was not explicitly considered, whereas the irregular component of AE_Budyko contained distinct variations in humid and sub-humid areas.


INTRODUCTION
Actual evapotranspiration (AE) is the process by which the Earth's surface water enters the atmosphere through evaporation or evapotranspiration [1]. Land AE plays an essential role in hydrology and climate change, as it links the water and energy cycles [1,2] and affects water-energy balances significantly. On average, AE accounts for about 60% of the Earth's total land rainfall [3]. Meanwhile, AE as latent heat flux is a central component of the heat transfer. A drop in evaporative cooling generally causes higher surface air temperature [4]. Global warming is expected to increase AE rates, further accelerating warming via water vapor feedback, since water vapor is a potent greenhouse gas [5]. It should be noted that, however, AE is influenced by both local and global factors [1,6]. The response of AE to global warming, if fully realized, can vary regionally across a continent [7].
AE is well known to be one of the most difficult variables to observe [8]. For now, only relatively few flux tower sites operate in China, and the period is concise [9]. The sparse observations are difficult to capture long-term variations in regional AE. Numerous theoretical and empirical methods have been proposed for estimating AE indirectly. For example [10], developed a data-oriented AE product by integrating a variety of data sources in a model tree ensemble approach, demonstrated that the majority of China (78%) showed an increasing trend in AE. [11], using a complementary relationship model, showed that annual AE increased significantly in western and northeastern parts of China. Recent improvements in retrospective analyses (or reanalyses) also enable AE estimates to be obtained at multiscale [12]. AE variability largely depends on local microclimatic conditions in a given region [13]. Various methods or models generally provide inconsistent AE estimates [11]. AE variability over China under global warming remains highly uncertain.
Due to the influence of the monsoon, the annual cycle is the major component of many climate variables in China [8]. Nevertheless, seasonality in time series can obscure movements of other parts that are operationally more important for climatological analysis. The annual cycle in most climate studies is assumed as a regular and recurring change with no interannual variation, whereas in reality climate change is reflected not only in the annual mean shift in climate variables, but also in their annual cycles [14]. Therefore, the seasonal, interannual, and interdecadal variations in AE and the associated drivers should be studied separately. Many deseasonalization methods and procedures have been proposed to deal with seasonal adjustment [14,15]. The widely used X-12-ARIMA method allows for interannual variations in terms of seasonality [15] and has the advantage of diagnosing AE variability.
Various AE estimation methods or models give inconsistent trends, a single AE dataset is insufficient [12,16]. Therefore, we examined AE variability over China derived from three sources, namely the Budyko framework, a complementary-relationshipbased product [11], and a weighted average of six reanalyses. The X-12-ARIMA method was then used to comprehensively evaluate the performance of these three datasets in reproducing the trend, seasonal, and irregular variations in AE for China's six climatic areas. The six areas are divided according to the aridity index, and coincide roughly with China's climatic characteristics. Furthermore, the differences in the contributing factors to AE variations were determined. Our results mainly focused on the recent changes in and attributions of regional AE.

Actual Evapotranspiration Data
Reanalyses (or retrospective analyses) have been produced by various institutes. In this study, AE data from six different reanalyses are used-namely, ERA-Interim reanalysis (ERA-I) from ECMWF [17], MERRA and MERRA2 from the NASA Goddard Space Flight Center [18], the Japanese 55-years Reanalysis (JRA-55), NCEP/NCAR reanalysis I (NCEP-R1), and NCEP/DOE reanalysis II (NCEP-R2) (available online https://rda.ucar.edu/). These reanalyses were generated by various forecast models, assimilation systems, or input datasets. Different reanalyses provide a measure of uncertainty in estimating the AE variability. For detailed descriptions of the selected reanalyses, refer to [12,13]; or the original dataset citations.

Meteorological Data and NDVI Data
To clarify the primary cause for variations in AE, we analyzed the relation between AE and the climate variables, which include surface net solar radiation derived from ERA-I, mean temperature, maximum and minimum temperatures, diurnal temperature range, cloud cover, potential evapotranspiration, precipitation, wet-day frequency, and vapor pressure derived from the University of East Anglia CRU TS4.00 version dataset (available online https://crudata.uea.ac.uk/cru/data/hrg/ cru_ts_4.00/; [20]. In addition, the 8 km Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index (NDVI) dataset was derived from the Advanced Very High Resolution Radiometer (AVHRR) sensor (available online https://www.ncei.noaa.gov/data/), and was used to estimate vegetation coverage [21]. Each dataset has a different time span, but all cover our analysis period of 1982-2015. For easy comparison and for computation purposes, all datasets were interpolated into a 0.5°× 0.5°grid using bilinear interpolation.

Methods
The X-12 Seasonal Adjustment Procedure Seasonal adjustment corresponds to an estimate of seasonal variation and its elimination from a time series [15,22]. Seasonally adjusted data can be used to better visualize the long-term development of the data series. We adopted a widely-used seasonal adjustment method, which is known as the X-12-ARIMA-based filter, to give a richer diagnosis and clearer interpretation of climate variations in AE. To be more precise, the original series (X t ) can be decomposed into three components, as where T t , S t , and I t are defined as the trend, seasonal, and irregular components, respectively [22]. The trend component (T t ) captures the direction of the time series, thus the nonparametric Mann-Kendall test was applied to detect its linear trend. Irregular fluctuations were generally due to unpredictable and unexpected factors. For more sophisticated analyses, it was preferable to include both the seasonally adjusted data and its seasonal component to obtain maximum information from the data [15].

Weighted Average of AE From Six Reanalyses
Reanalysis AE is generally estimated from bulk flux formulas with inputs of surface temperature, wind, and surface air temperature and humidity [16]. Discrepancies in AE among reanalyses can be substantial [12,13]. An intercomparison of different datasets is a key step in reducing uncertainties in AE estimates. After seasonal adjustment of AE, multiple linear regression, a sophisticated and quantitative analysis method, was used to analyze consistency across the six reanalyses.
For the X-12 trend, seasonal, and irregular components of AE from reanalysis X 1 , a multiple linear regression was performed with the remaining five reanalyses (X i , i 2, . . . , 5), respectively, that is, where β i (i 0, . . . , 5) are the parameters generally estimated by least squares. Based on this, the corresponding squared multi-correlation coefficient (R 2 i ) was used as the weight coefficient of each reanalysis. Then we calculated the weighted average of AE (AE_WAR) from the six reanalyses using Eq. 3, Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 718771 AE WAR where the subscript i of R 2 i represents different reanalysis data mentioned above.

AE Estimation Based on the Budyko Framework
In the Budyko framework, long-term AE (AE_Budyko) is largely controlled by the ratio of energy supply (potential evapotranspiration, PET) to water supply (precipitation, PRE). Several previous studies have developed various but somewhat similar coupled water-energy balance equations on the basis of the Budyko hypothesis [23]. Because of its strong theoretical basis and low data requirements, the Budyko equation has been widely used [24]. As pointed out by [25]; various types of Budyko-like equations have given broadly similar results. One of the most popular Budyko equation is expressed as [21]: where n is a catchment-specific parameter that modifies the Budyko curve [24,25]. China occupies a vast land territory, thus the parameter n shows a large spatial variation (ranging from 0.4 to 3.8) [25]. It should be noted that higher values of n denote a higher estimate of AE for a given PRE and PET in general [24]. To simplify the Budyko framework, the parameter n is traditionally set to be 1.8 [13,23]. In the absence of better knowledge, we assume here that n remains constant over time.

Climatic Characteristics of China
Climate variations in China are complicated because of the influence of various factors. Precipitation varied widely across space ( Figure 1D). Generally, wetter environments had a higher vegetation cover ( Figure 1C), indicating that long-term vegetation cover was sensitive to water availability. Considering geographical complexity, the response of regional AE to a warming climate should be investigated separately [13].
Meanwhile, the quality of the AE products may differ spatially as well. We thus differentiated the climatic areas in China and performed an attribution analysis on AE for each area. According to the aridity index (AI) [26], the Budyko hypothesis generally separates catchments into energy-limited (AI < 1) and water-limited (AI > 1) areas [27]. As shown in Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 718771 Figure 1A, AI had a similar spatial distribution to precipitation ( Figure 1D) and vegetation cover ( Figure 1C). Consequently, the definition of AI was deemed reasonable and reliable. According to this index, China could be roughly divided into six climatic areas, namely, the hyper-arid (H-arid), arid, semi-arid (S-arid), dry subhumid (DS-humid), sub-humid (S-humid), and humid areas [13]. Figure 1B shows the spatial distribution of climatological mean AE during 1982-2015, calculated from the average of six reanalyses. Generally, AE across China decreased from south to north and from east to west. For the mean annual water balance, the moisture supplied by precipitation was sufficient in most regions. In contrast, AE could be higher than precipitation in northwestern China, where glaciers have a direct impact on the water resource system [39]. Precipitation is used to measure water availability in the Budyko framework. It means that AE_Budyko may be biased low in some areas (discussed below). Figure 2 shows the seasonal component of AE for the six climatic areas. AE showed strong seasonality with different amplitudes over time. In comparison, seasonal variations in AE revealed by different datasets were generally consistent for most areas. The largest discrepancies were observed for the hyper-arid area ( Figure 2A). Overall, AE_WAR agreed better with AE_CR (correlation 0.926) than with AE_Budyko (correlation 0.856). Modest differences also existed in the seasonal component for the arid area ( Figure 2B). Seasonal variation is an essential aspect of climate change. As shown in Figure 2, seasonal components of AE in the hyper-arid and arid areas exhibited substantial interannual variations, which is a nonnegligible part of the original time series. It suggested that the conventional approach assuming fixed seasonality did not apply in these areas.

Multi-Timescale Variability of AE
The trend components (T t ) of AE are shown in Figure 3. Compared to AE_CR, AE_WAR agrees better with AE_Budyko, as evidenced by the correlation coefficients shown in Figure 4 (0.589, 0.623, 0.708, 0.540, 0.387, and 0.288 for the six areas, respectively). Despite the general similarity, disparities existed among the three datasets, partly owing to their differing linear trends ( Figure 5). The statistical significance of the Mann-Kendall test for the increasing trend in AE_WAR exceeded the 99% confidence level in all areas ( Figure 5). This was consistent with the traditional view that AE is likely to increase under global warming [7]. AE_Budyko also showed an increasing trend in most areas, whereas a significant negative trend was found for the hyper-arid area and no Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 718771 significant trend existed for the arid area ( Figure 5). AE_CR displayed a significantly decreasing trend in all areas, except for the arid area. The spatial patterns of the linear AE trends, and indeed the signs of the trends, differed between the datasets. Compared to AE_WAR and AE_Budyko, AE_CR showed a significant discrepancy in the trend component ( Figure 4). Specifically, AE_CR for many areas peaked in 1997 (Figure 3), increasing during 1982-1997 and declining after 1997. We suspected that the aforementioned weak correlation might have been caused by the effect of the spurious changes in AE_CR around 1997. A correlation analysis between AE_CR and AE_ WAR was performed for 1982-1995 and 1999-2015, respectively. Correlation coefficients were highly improved during 1982-1995. The correlations between AE_CR and AE_WAR were 0.500, 0.692, 0.702, 0.712, 0.519, and 0.413 (p < 0.001) for the six regions, respectively. The relationship between AE_CR and AE_Budyko showed roughly similar behavior. A comparison of various datasets revealed substantial uncertainties in the trend component of AE_CR, which was likely overestimated around 1997. It is therefore questionable whether a calibration-free AE estimation method without the use of precipitation or runoff rates could predict long-term changes in AE.
AE_WAR increased from the mid-1990s to the early 2000s and more recently in the hyper-arid area ( Figure 3A); this was generally consistent with the AE_Budyko trend. Similar results were found for the arid area ( Figure 3B), where AE_WAR also had an increasing trend in the late 1980s. At a regional scale, AE_WAR agreed better with AE_Budyko for the semi-arid area ( Figure 3C), with both indicating a distinct increase in AE from 1982 to the early 1990s. During the same period, an increasing trend in AE_WAR was observed for the dry sub-humid area, whereas AE_Budyko showed no discernible trend ( Figure 3D). Another significant shift in both AE_WAR and AE_Budyko occurred from the late 2010s onwards. Considerable differences in the trend components of AE among the three datasets were found in the sub-humid ( Figure 3E) and humid ( Figure 3F) areas. For these two areas, AE_WAR had a significantly increasing trend along with substantial interannual variability, which agreed with AE_Budyko. As noted above, both AE_WAR and AE_Budyko agreed well with AE_CR during 1982-1996. However, variations in AE obtained from the three datasets were not consistent and exhibited conflicting trends during 1997-2015. At a regional level, multiple factors determined how much water evaporated from the land [1]. The irregular component of AE showed high variability ( Figure 6). Considerable inconsistencies existed between the three datasets. Specifically, AE_WAR showed reasonable agreement with AE_Budyko for the drylands (i.e., hyper-arid, arid, and semiarid areas). By contrast, only a weak negative correlation existed between AE_WAR and AE_Budyko for the humid area ( Figure  5F). One possible reason for this was that changes in both rainfall and soil water storage significantly influenced the interannual variability of the hydrological responses [28]. For example, vegetation and the interactions between climate seasonality and soil water storage changes have also been  found to play important roles [25]. The magnitude of any change in the interannual variability of AE caused by variations in catchment vegetation was related to pre-and post-change vegetation types. Compared to the drylands, the vegetation cover in wetter environments tended to be denser ( Figure 1C). Because of the mutual interactions among climate, hydrological processes, and land surface characteristics, AE had a nonlinear relationship with precipitation and potential evapotranspiration [24,29], which led to deviations in the Budyko relationship [28]. As a result, the performance of AE_Budyko in reproducing the irregular component of AE remains doubtful for the humid area ( Figure 6F).
Agreement between AE_WAR and AE_CR was observed for the humid and sub-humid areas ( Figure 4). As the above mentioned, a calibration-free nonlinear complementary relationship model was utilized for monthly AE_CR taking air and dew-point temperature, wind speed, and net radiation as inputs [11,30]. AE from surfaces with abundant moisture was mainly controlled by energy conditions [27,31]. Atmospheric evaporative demand was primarily driven by two major components, namely radiative and aerodynamic components, both of which were associated with the AE_CR inputs. Therefore, the irregular component of AE_CR performed better than that of AE_Budyko for the humid and sub-humid areas (Figure 4).

Causes and Implications of Changes in AE
The response of AE to climate variables varies spatially (Figure 7). The causes of variations in AE are closely related to local climatic characteristics. The trend component of AE showed a strong positive correlation with precipitation, wet-day frequency, and vapor pressure (p < 0.001) (Figure 7) in the hyper-arid, arid, and semi-arid areas. This was consistent with the convention that AE is mainly controlled by water conditions in drylands. In waterlimited environments, a substantial fraction of the variations in precipitation become variations in AE. Although precipitation was the dominant factor for the trend component of AE in the drylands, the increase in precipitation (Figure 8) was lower than that in AE (Figure 4) according to the Mann-Kendall test.
In fact, increased wet-day frequency could promote AE increase, even if there is no change in total precipitation. Besides, changes in runoff and water storage should be considered in the amount of moisture available for AE [31,32], especially in drylands. Human activities and climate change have intensely influenced the ecohydrological pattern of many basins. For example, the Tarim River Basin is the most heavily glacierized watershed in arid northwestern China. The water supply from the mountains to the Tarim River Basin is favorable in recent decades; however, local human activities (such as irrigation and domestic water use) have reduced its runoff [33]. [34] also indicated that the influence of irrigation on AE variability became increasingly evident with the increase of irrigation in the Tarim River Basin. Accordingly, the increased water supply from catchment storage likely resulted in increased AE in the drylands. AE_WAR and AE_Budyko were significantly correlated (p < 0.001) with diurnal temperature range, net solar radiation, and cloud cover in the drylands. The Mann-Kendall test indicated that diurnal temperature range decreased significantly in the drylands, mostly because of a more significant increase in minimum temperature relative to maximum temperature ( Figure 8). This decrease in diurnal temperature range was mainly attributed to increased cloud cover, precipitation, and soil moisture [35]. Diurnal temperature range can represent the combined influence of these variables and hence has a strong relationship with AE [36].
Increased AE can modulate the climate by reducing the sensible heat flux, enhancing the air humidity, increasing the minimum temperature, and decreasing the maximum temperature [35]. It means that increased AE is a possible reason for the increase in minimum temperature has been higher than that in maximum temperature, that is, the decrease in diurnal temperature range might have resulted from the increased AE (Figure 8). The correlations ranged more widely between the datasets in wetter areas. Vapor pressure was an essential factor influencing AE in the dry sub-humid area ( Figure 7D). Increased AE led to increasing vapor pressure, thus increased vapor pressure could be interpreted as a sign of increasing AE [6]. Similar relationships were found in the sub-humid area ( Figure 7E). Overall, the simultaneous increases in precipitation and wet-day frequency increased AE in the dry sub-humid area ( Figure 7D). Increased minimum temperature was likely an essential cause of the increased AE in the sub-humid area ( Figure 7E). For the humid area, the trend component of AE was less sensitive to precipitation than that in the drylands. AE increased under the integrated influence of potential evapotranspiration, mean temperature, and maximum and minimum temperatures ( Figure 7F). Changes in AE were generally controlled by energy supply rather than precipitation.
In general, AE was driven mainly by climatic factors, regulated by land surface condition, and limited by water supply [19]. The trend component of AE_CR showed little correlation with climatic factors. An explanation is that AE_CR may fail in reproducing the long-term variation in AE. In contrast, AE_Budyko showed reasonable consistency with AE_WAR, demonstrating that AE_Budyko is more reliable in terms of long-term, large-scale changes in AE. Figure 9 illustrates the influencing factors of the irregular component of AE at regional scales. As expected, precipitation had the strongest influence on AE in the hyper-arid, arid, and semi-arid areas. In the drylands, precipitation is temporarily stored in soil and then totally evaporates into the atmosphere [1]. In contrast, AE_WAR and AE_CR were significantly negatively correlated with precipitation in the humid area ( Figure 9F), where AE was more sensitive to energy supply than water supply ( Figure 1A). The irregular component of precipitation made a much smaller contribution to AE. Furthermore, precipitation usually coincides with an increase in cloud cover and a decrease in radiation, resulting in a reduction of energy supply to AE ( Figure 9F). For this reason, precipitation had a negative effect on AE in the humid area. Similarly, positive correlations between AE and cloud cover dominated in drylands, but significantly negative associations existed between AE and surface net solar radiation.
AE was significantly negatively correlated with potential evapotranspiration in the drylands (Figure 9). [29] suggested that a complementary relationship between AE and potential evapotranspiration occurs in non-humid areas because these two factors are correlated via precipitation. In contrast, AE and potential evapotranspiration had a strong positive correlation in wet areas ( Figures 9E,F). The inclusion of the regional dimension of AE drivers allows the scenarios both of increasing AE and potential evapotranspiration in water-rich areas and of increasing AE and decreasing potential evapotranspiration to be encompassed in the drylands (Figure 9).  In the humid area, surface net solar radiation was the most crucial variable regulating the irregular component of AE_WAR ( Figure 9F), which was reasonable, as radiation is the energy source for AE. Overall, the result of AE_WAR matched that of AE_CR ( Figure 9F). [37] also reported that in humid climates, such as northeast India, solar radiation supplied most of the energy required for water to change from a liquid to a vapor. Furthermore, AE was strongly correlated with diurnal temperature range. [1] pointed out that the variability in diurnal temperature was consistent with that in surface incident solar radiation at monthly to decadal timescales. Since surface incident solar radiation heats daytime air only, its variation was related to changes in diurnal temperature.
AE is a complex process that regulates land-atmosphere interactions. Increased water supply from catchment storage could be a primary contributor to increased AE. Nevertheless, accurate detection and attribution of complex variability in AE at multi-time scales remain challenging [33,38]. The implications of variations in AE remain unclear.

CONCLUSION
According to the aridity index, China was divided into six climatic areas. This study analyzed regional AE variability across China using data from the weighted average of six reanalyses (AE_WAR), the Budyko framework (AE_Budyko), and a complementary-relationship-based AE product (AE_CR). Since seasonal variation complicates the detection of changes in AE, the X-12-ARIMA method was used to study the trend, seasonal, and irregular components of AE as well as the significant drivers.
According to the trend component, AE increased from the mid-1990s to the early 2000s and more recently in the hyperarid and arid areas. An increasing trend in AE was also observed from 1982 to the early 1990s in the semi-arid and dry sub-humid areas. In the sub-humid and humid areas, AE had a significantly increasing trend along with substantial interannual variability. Increased precipitation and water supply from terrestrial water storage contributed significantly to increasing AE in drylands. The simultaneous occurrences of increased precipitation and wet-day frequency caused high AE in the dry sub-humid area. In contrast, increased AE in the sub-humid and humid areas could be explained by increased energy supply and precipitation.
Precipitation had the strongest influence on the irregular component of AE in the drylands. In contrast, AE_WAR and AE_CR were significantly negatively correlated with precipitation in the humid area. Negative correlations between AE and potential evapotranspiration occurred in the drylands, whereas positive relationships existed in the sub-humid and humid areas. The irregular component of AE was positively correlated with vapor pressure across China. Positive relationships between AE and cloud cover dominated in drylands, but significantly negative associations existed between AE and surface net solar radiation.
Considering data availability, the seasonal variations in AE revealed by the different datasets were generally consistent. The largest discrepancies occurred in the hyper-arid area. Compared to AE_WAR and AE_Budyko, the trend component of AE_CR showed inconsistent behavior across many areas, which was mainly caused by spurious changes around 1997. AE_CR did not consider soil moisture explicitly, which may have limited its ability to detect AE trends. The Budyko model is useful in largescale hydrological research; however, the irregular component of AE_Budyko remains dubious for humid and sub-humid areas because of the complex interactions between climate, vegetation, and soil.
Water loss through AE is a major consideration in the design and management of water supply reservoirs. Efforts to evaluate AE in natural settings are made difficult by spatial heterogeneity in soil and unevenly distributed vegetation in addition to other biophysical processes. This study thoroughly compared different AE products for various climatic areas of China. Such knowledge is meaningful to hydrological, climatological, and ecological research.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author. The GIMMS NDVI dataset, the Climate Research Unit datasets, and the six reanalyses used in this study can be freely accessed from the following websites: https://www.ncei.noaa.gov/data/, https:// crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.00/, and https://rda. ucar.edu/, respectively.

AUTHOR CONTRIBUTIONS
TS and GF conceived, and carried out the research, led the data analysis, and wrote the manuscript. TF and BH conceived the research, designed the analysis, and provided comments on the manuscript. ZH, ZQ, and WH provided comments on the manuscript.